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Abstract 

A universal particle velocity based algorithm for simulating hydraulic fractures, previously pro¬ 
posed for Newtonian fluids, is extended to the class of shear-thinning fluids. The scheme is not 
limited to any particular elasticity operator or crack propagation regime. The computations are 
based on two dependent variables: the crack opening and the reduced particle velocity. The appli¬ 
cation of the latter facilitates utilization of the local condition of Stefan type (speed equation) to 
trace the fracture front. The condition is given in a general explicit form which relates the crack 
propagation speed (and the crack length) to the solution tip asymptotics. The utilization of a mod¬ 
ular structure, and the adaptive character of its basic blocks, result in a flexible numerical scheme. 
The computational accuracy of the proposed algorithm is validated against a number of analytical 
benchmark solutions. 


1 Introduction 

The phenomenon of a hydraulically propelled fracture that propagates in a solid material is present in 
many natural processes and in technology [3 031GH ED- The multiphysical character of the problem 
makes mathematical and numerical modeling of hydraulic fractures (HFs) an extremely challenging 
task. The main computational difficulties stem from: i) strong non-linearities resulting from interaction 
between the solid and fluid phases, ii) singularities in the physical fields, iii) moving boundaries, iv) 
degeneration of the governing equations at the singular points of the domain, v) leak-off to the rock 
formation, vi) pronounced multiscale effects, vii) complex geometry, and others. 

The first simplified mathematical models of hydraulic fractures were introduced in the 1940s and 
1950s [ffil [24j [26j (28] [59]. These works, together with some later studies, evolved into basic ID models: 
i) the PKN model EDI iS], ii) the KGD model [19, 2Tj, iii) the radial or penny shaped model [551 . 
These models were used for decades in practical applications to design the hydrofracturing treatments. 
However, a continuous increase in both the number and size of fracking installations necessitated further 
development of the mathematical models describing the process, which resulted in the introduction of 
the so-called pseudo 3D models [IS] , planar 3D models (4] [9J [65: and recent attempts to develop fully 
3D models Eunnsacn]. A comprehensive review of the topic can be found in [2]. 
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Simultaneously, the studies of the basic solution features pertaining to the underlying physics of the 
process have been carried out. Special attention has been devoted to the near-tip region of the fracture. 
The asymptotic behaviour of the solution for various hydraulic fracture models and crack propagation 
regimes has been identified in a number of works - see e.g. mmmm and references therein. 
The importance of the multiscale character of the problem has been fully recognized. It is now well 
understood that the global behaviour of the hydraulic fracture is primarily controlled by its near-tip 
region [18] . As a result of the moving boundary of the domain, the governing equations and boundary 
conditions degenerate at the crack tip, which makes tracing the fracture front an extremely difficult task 
[53] . Moreover, the solution tip asymptotics essentially depend on the fracture propagation regime. Here, 
four basic modes have been identified H8IEZ]: i) leak-off dominated, ii) storage dominated, iii) toughness 
dominated, iv) viscosity dominated. All these factors clearly indicate that a proper understanding of the 
solution structure and the correct application of the tip asymptotics in the computational schemes are 
crucial for the accuracy, credibility and robustness of computations. In the recent paper by Lecampion et 
al. [38] it was shown that the algorithms which utilize the appropriate multiscale hydraulic fracture tip 
asymptote perform much better than those that do not apply it. The advantages of employing enhanced 
asymptotic representation in the numerical algorithm have been demonstrated in [17] 12711301 (special tip 
element), [2] [5] EH 1221[23]133 EZ] (BEM and XFEM formulations). 

Another crucial element affecting the credibility and efficiency when numerically modeling hydraulic 
fractures is the mechanism of fracture front tracing. As in any other problem involving a moving 
boundary, a proper Stefan type condition should be imposed. However, its direct utilization in HF is 
a numerically challenging task M- Alternatively, in the case of ID problems (PKN, KGD and radial 
models), the fluid balance equation can be used to trace the crack tip position. The non-local nature 
of such a condition has a clear advantage from the computational point of view, in comparison with 
any local counterpart. When properly employed in the numerical scheme, it is equivalent to the Stefan 
condition and allows one to obtain accurate results for a reasonable computational cost. As for the 
Stefan type boundary condition, its applicability and importance in hydraulic fracture problems was 
pointed out in [25]. Two decades later its relevance was rediscovered in [35], where the condition itself 
was called the speed equation. 

The situation is quite different when one analyzes 2D or 3D models of hydraulic fractures. Thus, it 
is of vital importance to utilize a method based on the local condition to establish the fracture front 
location. To this end an implicit level set method was developed in [52]. The technique was successfully 
implemented in [53]. Recently, an explicit level set scheme was proposed in [44] . It combines the 
information in the solution tip asymptotics with the speed equation. 

Such an approach for the ID PKN model has already been presented in [32] 1461 : (j 51 while its gen¬ 
eralization, introducing explicit relations between crack tip asymptotics and crack propagation speed 
for all classical ID models, can be found in [651 . It was shown that, when properly combined with 
other numerical devices, the algorithm represents an accurate and universal tool, which allows one to 
circumvent the main problems found in classical approaches used to establish the crack tip position. In 
this paper, we will consequently use the methodology developed in m- 

Appropriate modeling of the interaction between solid and fluid phases necessitates a thorough in¬ 
vestigation into the underlying physical mechanisms of fluid flow inside the fracture. One of the most 
important problems here is a proper description of, and accounting for, the rheological properties of the 
fracturing fluid. Depending on the environment in which the fracking process is carried out, various types 
of fluids are in use. In particular, in low permeability reservoirs (small leak-off to the rock formation) one 
usually employs water with some additives. In such a case the fluid can be modeled as a Newtonian one. 
Most of the fracturing fluids utilized in the conventional reservoirs exhibit non-Newtonian properties 
that can be described by a power law [6J [T5J (33, 3]. 135] j6J[ : 

r = My”. (1) 
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Here r denotes the shear stress, 7 stands for the shear strain rate, n is the ’’fluid behaviour index” and 
M the ’’consistency index”. The range of 0 < n < 1 refers to the so-called shear-thinning fluids, while 
n > 1 pertains to the shear-thickening case. For n = 1, one has the Newtonian fluid with the consistency 
index equal to the dynamic viscosity: M = p. A special case when n = 0 pertains to the perfectly 
plastic fluid, where M = ay is the yield stress. 

However, as pointed out in EH, such a simple two-parameter model cannot correctly describe the 
real fluid behavior at low and high shear rates. In such cases more advanced approximations, like the 
Carreau model dim are appropriate. Unfortunately, a closed-form solution for the Poiseuille type 
flow in a channel with parallel walls cannot be obtained with the Carreau law. Thus, the lubrication 
theory is no longer in use here. In order to circumvent this deficiency, the so-called truncated power 
law was recently proposed in [36] . This four-parameter approximation uses cut-off viscosities at low and 
high shear rates, with the power law dependence in the intermediate range. It provides a satisfactory 
imitation of the Carreau model and simultaneously allows one to obtain the closed-form solution for the 
flow between parallel plates. A three-parameter rheological model for the Hershel-Buckley fluid type 
[25] was used in [ 22 ], where the power law 0 was combined with a lower cut-off when the shear stress 
becomes equal to the value of the yield stress. A closed-form solution for the Poiseuille type flow in the 
channel was embedded in the PKN model with the Carter leak-off law. 

Nevertheless, as far as the simulation of hydraulic fractures for non-Newtonian fluids is concerned, 
most of the results available in literature refer to the pure power law dependence of the form ([Tj). 

An analytical solution of the PKN model for shear-thinning fluids was obtained in m- The author 
reformulated the problem in terms of new dependent variables, using the power of the crack opening (a 
modified opening) and the fluid velocity. A self-similar solution in the form of an infinite power series was 
given for the case of a predefined leak-off. Particular solutions for an impermeable rock (zero leak-off) 
for two limiting values of the fluid behaviour index, n = 0 and n = 1 (perfectly plastic and Newtonian 
fluid), were later used in [43] to draw conclusions concerning the equivalence of different shear-thinning 
fluids with respect to their hydrofracturing action. In mi, the authors analyzed the PKN problem for 
the zero leak-off case. They introduced two simple analytical solutions: for constant fluid velocity along 
the fracture length and for constant crack volume in time. A numerical algorithm based on reducing 
the problem to a dynamic system of the first order was proposed. A number of numerical results for 
shear-thinning and shear-thickening fluids were presented. 

A penny-shaped fracture induced by a power-law fluid was analyzed in [3]. The governing equations 
of the problem were derived and explicit time-dependent particular solutions were proposed for both 
Newtonian and shear-thinning fluids. Parametric sensitivity studies were conducted to estimate the 
influence of both the fluid viscosity and the leak-off (Carter law) on the fracture evolution. 

The KGD model in viscosity dominated regime for impermeable rock was investigated in [Tj. The 
authors reduced the problem to a self-similar formulation and delivered semi-analytical solutions, in 
the form of power series, for both the Newtonian and a number of shear-thinning fluids. In [30; the 
KGD model for a viscosity dominated regime was analyzed in the case of permeable rock. A finite 
element based algorithm with no explicit fracture front tracing mechanism was proposed to find the 
solution. The credibility of computations was verified against the reference solutions from [T, . The KGD 
model for a finite material toughness was investigated in m ■ The author proved that for shear-thinning 
fluids the crack propagation regime evolves in time from toughness dominated to viscosity dominated. 
Corresponding approximate self-similar solutions for these limiting cases were delivered in the form of 
power series. A numerical algorithm for the transient solution was also proposed. 

In this paper we analyze the problem of a ID hydraulic fracture driven by a non-Newtonian power 
law fluid whose rheology is described by equation ([!]). Two classical hydrofracturing models, the PKN 
and the KGD (for the viscosity dominated and toughness dominated regimes), are under consideration. 
A unified numerical algorithm, being a direct extension of the one introduced for Newtonian fluids in 
[69] . is proposed. The basic assumptions of the numerical scheme are: i) a proper choice of independent 
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and dependent variables (reduced particle velocity), ii) a mechanism of fracture front tracing based on 
the Stefan condition, utilizing a universal explicit relation between the crack propagation speed and the 
coefficients in the asymptotic expansion of the crack opening, iii) proper regularization techniques, iv) 
an improved temporal approximation, v) modular algorithm architecture. 

Some of the aforementioned elements were already presented in [32, j|f5, JJS], where one can find a 
broad analysis on the application of proper independent and dependent variables and regularization tech¬ 
niques. The theoretical background of the mechanism of fracture front tracing together with a detailed 
description of the techniques of its application was given in [59] . Adaptation of the numerical scheme 
from [69] to the class of shear-thinning fluids was primarily related to the different form of the Poiseuille 
equation and necessitated: i) updating the asymptotic expansions, and the interrelations between them, 
for respective dependent variables, ii) modifying the regularization techniques, iii) development of new 
analytical benchmark solutions utilized to verify the algorithm. The accuracy and efficiency of compu¬ 
tations, investigated against dedicated analytical benchmarks, demonstrate the ultimate performance of 
the proposed numerical scheme. 

The structure of the paper is as follows. In Section [2j a general formulation of the problem is given. 
After subsequent normalization the complete system of equations is reformulated in terms of a new pair 
of dependent variables: the crack opening and the reduced particle velocity. The mechanism of fracture 
front tracing is introduced. The employed methodology is based on proper utilization of the solution 
tip asymptotics. Section [3] is devoted to the time-independent (self-similar) problem formulation. The 
self-similar variant of the universal algorithm is introduced and thoroughly tested for the considered 
hydraulic fracture models. The performance of the algorithm is verified against analytical benchmarks 
presented in Appendix A, as well as reference solutions found in literature. A new approximation for 
the classical problem of the KGD model in viscosity dominated regime is delivered. In Section [4] the 
transient version of the algorithm is presented. Again, its performance is validated in a number of tests. 
The main peculiarities of the scheme are shown. The remaining part of the paper contains discussion of 
the results and final conclusions. Additional material is given in the appendices. 


2 Problem formulation 


Let us consider a rectilinear symmetrical crack of length 21 situated in the plane x £ [—1,1], The crack 
is fully filled by a non-Newtonian liquid injected at the middle point (x = 0) at a known rate qo(t). As 
a result, the crack front (x = ± l ) moves and the crack length, l = l[t ), is a function of time. Below we 
present a classic set of governing equations for the ID formulation of the problem, which can be found 
in e.g. [mm eh for various hydraulic fracture (HF) models. As usual, due to symmetry of the problem, 
we restrict our analysis to one of the symmetrical parts of the crack x £ [0, l(t)]. 

The mass conservation principle is expressed by the continuity equation: 

— + ^-+®=0, t > t 0 , 0 < x < l{t), (2) 


while fluid flow inside the fracture is described by the Poiseuille equation: 


q 


n 


l 

M 7 


2n+l 


dp 

dx 


(3) 


In the above equations w = w(t,x) is the crack opening, q = q(t,x ) stands for the fluid flow rate, 
p = p(t,x ) (p = pf — (Toj - confining stress) refers to the net fluid pressure. The constant M' in 
the Poiseuille equation is a modified fluid consistency index M' = 2” +1 (2n + 1 ) n /n n M. The function 
qi = qi(t, x) on the right hand side of ([2] is the volumetric fluid loss to the rock formation in the direction 
perpendicular to the crack surfaces per unit length of the fracture. In general, this function is itself a 
part of the solution, for which an additional diffusion problem in the surrounding rock needs to be 
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solved. Such a task goes beyond the scope of our analysis and thus we assume that the leak-off function 
is predefined. 

Fluid equations ([2]) - ([3]) are supplemented by the equation describing rock deformation under applied 
hydraulic pressure. Here we use the following symbolic relation: 


p(t,x) = Aw(t,x), 0 <x<l(t), 


( 4 ) 


where the operator A refers to the chosen HF model (related to the fracture geometry). In this paper 
we consider two hydraulic fracture models: 


• the PKN model |JT] 




>4iu> = k\w, 

( 5 ) 

• the KGD model jU QT 




k 2 f l(t) dw{t, s)s 

AlW =2 nl ds 

(6) 


For the classical PKN model, the constant k\ can be found in m for an elliptical crack of height h, 
where E and v are the elasticity modulus and Poisson’s ratio respectively. The multiplier k 2 in the KGD 
model follows, for example, from [1§1 IHT)] : 


2 E 

7r/l(l — u) 1 


k l = -T7-, - k 2 = 


E 

(i^y- 


The inverse operators for © © 


are: 


A p= i-p, 
k i 


4 f 1 ^ 

A 2 p = —r- / pit, s) In 
Jo 


^/P(t)-x 2 + ajl 2 {t)-s 2 


y/l 2 (t)~x 2 - ^P{t)-s 2 


ds. 


The form of the operators © and ([9]) implies that: 


dw 

dx 


(t, 0) = 0. 


( 7 ) 

( 8 ) 
( 9 ) 


( 10 ) 


The foregoing equations are supplemented by the influx boundary condition at the crack mouth 


(x 


0 ): 


q{t,0) = q 0 (t), 


( 11 ) 


and two boundary conditions at the crack tip: 

w(t,l(t)) = 0, q(t,l(t)) = 0. (12) 

In this paper we assume non-zero initial conditions: 


Z(0) = l 0 , w(0,x) = w 9 (x), xG(0,l o ). 


( 13 ) 


Note that (131 in its general form accounts for a preexisting hydraulic fracture. 

In the toughness dominated KGD model the fracture evolution is analyzed in the framework of Linear 
Elastic Fracture Mechanics (LEFM), and the following propagation condition holds: 


Kr(t) = K IC , 


(14) 
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where Kjc is the material toughness [55] , and Kj is the stress intensity factor (SIF) defined as 'TTfl : 


Ki _o./W f' (t) p(t,s)ds 


71 Jo \/l 2 (t) - s 2 


As a result the crack opening asymptote yields 


(t,x) ~ 8 


(15) 


(16) 


Finally, the global fluid balance equation, called also the solvability condition [2], can be written as: 
pZ(i) pt /*Z(t) rt 


r L \ l ) r* rw r* 

/ [w(t, x) — w 0 (x)\dx — / qo(t)dt+ / / qi(t, x)dtdx = 0. 

Jo Jo Jo Jo 


(17) 


The above equations constitute the classical model of hydraulic fracture lanang. 

In our analysis we shall use another dependent variable - the particle velocity, v. It describes the 
average speed of fluid flow through the fracture cross section and can be defined as: 


v(t, x) = 


q(t-,x) 
w(t , x) 


1 

= I- w 

1 M' 


n+1 


dp 

dx 


l/r 


(18) 


The advantages of applying this variable were shown in 169] . Moreover, its introduction is related 
to the mechanism of fracture front tracing utilized throughout this paper. Namely, provided that there 
is neither lag between fluid front and crack tip nor an invasive zone ahead of the fracture (the fluid 
front coincides with the fracture tip), the crack propagation speed equals the value of particle velocity 
at the crack tip (see discussion in nans esi). This allows us to formulate the relation for fracture front 
tracing as an evolution of the Stefan condition 123 EH, called also the speed equation 




1/n 


c=l(t ) 


(19) 


The technical implementation of condition (19) and advantages of such an approach are extensively 
discussed in [5^1. 


Remark 1. In general, (19) should be modified when the leak-off is described by the Carter law na. 


2.1 Normalization 


Let us normalize the problem by introducing the following dimensionless variables: 


i(ty 

q(t, x) = 


t 

t — yjn ’ 

w(t , x) ■ 

w(t , x) 

/ ’ 

C/* 

m = y, 

y/ n 

qi(t,x) = -T— qi(t,x), 


1/n 


1 , 


t.y™ , N m' 

(20) 

~pT<l(t,x) 

, p(t, X, 

) = — p(t, X 
K"m 

), v(t,x) = 

——v(t,x), t r = —, 
l* rCrn 


where x £ [0,1], n is the fluid behaviour index defined above ([Tj) , and we have that m = 1 — r, with r 
being the homogeneity order of the operator A m - The parameter m has value 1 or 2 for the PKN and 
KGD model respectively. This type of normalization has already been used in 1691 . It is not attributed 
to any particular influx regime, elasticity operator or asymptotic behaviour of the solution. 
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Using the transformations (201 and expressing the fluid flow rate in terms of the crack aperture and 
the particle velocity, one can rewrite the continuity equation |2]) as: 


dw L' _ dw 1 dtwv) 

~dt^Y x ~di + L~dT + qi = {) - 


The normalized particle velocity (18) yields: 


v = — = ——w 


-n+l 


dp 

dx 


1/n 


The notation for the elasticity equation Q is: 


with operators now defined as: 


A 2 W = 


p =Aw : 

Aiw = ui, 

1 1 r 1 dw(t,p) p 

2 tt L(t) Jo dr] x 2 - p‘ 


-dip 


while their inverses are: 


w = A 1 1 p = p 


= A 2 1 p = ~L(t) / p(t,rj) In 
7T Jn 


y/1 - X 2 + y/l - ' 
y/l — X 2 — \/1 — 


dp. 


Following [5^, we shall use an alternative form of (271 defined as: 

A^p = L{t ) [ K(rj, x )dp + -^=JL^K^l - x 2 , 

it J 0 dp V 7T v 


w = . 


where the kernel K(p : x) is: 

K(p , x) = (p — x) In 


yjl — x 2 + y/l — p 2 

y/1 — X 2 — y/l — p 2 


■ x In 


l + px + \/l — x 2 ^/l — p 2 
1 + px — \/l — X 2 \/l — p 2 


( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 


Kj = Kj c (l — v 2 )/(Ey/TJ) stands for the dimensionless toughness. Note that the fracture propagation 
condition (14) is already included in the formulation. Consequently, the asymptotic estimate (161 can 

(30) 


be rewritten as: 


\2tt v 


x —> 1. 


From definition (15) one can determine the dimensionless toughness as: 


K, = 2\ 


I L(t ) f 1 p(t, s)ds 


VT^' 


(31) 


The normalized global fluid balance equation gives: 


f [L(i)w(t, x) — L(0)w*(x)] dx — f qo{t)dt + f L(t) f qi(t,x)dxdt = 0, (32) 

Jo Jo Jo Jo 
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while the speed equation (19) converts to: 


Vo(t) 


dL 

dt 


1 

L 


w 


n+1 


dp 

dx 


- 1/n 


< OO. 


(33) 


Here, the value of vq( t) is bounded and equal to the crack propagation speed: vq( t) = v(t, 1). This 
assumption is valid if the crack opening and the leak-off functions behave asymptotically as: 


= w 0 (t)(l - x) ao + o({l-x) a °), qt =qio(t)(l-x) ri + o((l-x) TI ), 


1 , 


(34) 


where 1 + rj > a 0 (see discussion in [85]). Note also that relation (32) assumes that the process of the 
fracture growth is monotonic (£ 0 (f) = L'(t ) > 0). 

As previously, for the KGD model, the fracture symmetry condition holds: 


dw 

dx 


(f, 0) = 0 . 


The boundary conditions (11) |12j) assume the form: 

w(t,0)v(i,0) = qo(t), w(t, 1 ) = 0 . 


(35) 


(36) 


Note that the equivalent of condition @ 2 , when represented in terms of w and v, and accounting for 
condition (331, is satisfied automatically. 


The initial conditions are now as follows: 


m = ‘f, 


u)( 0 , x) = w*(x) = —w 0 (l 0 x), x £ [ 0 , 1 ]. 


(37) 


For simplicity, from now on we omit the symbol. 


2.2 Tip asymptotics and crack propagation speed 

Recently, a universal algorithm for numerical simulation of hydraulic fractures has been proposed [65 .. 
One of its main elements is the mechanism of fracture front tracing based on the Stefan condition. It 
also extensively utilizes the crack tip asymptotics. In this paper we are going to extend this approach to 
the case of shear-thinning fluids. To this end, below we provide the basic information on the asymptotic 
behaviour of the solution at the fracture front. 

As was shown in ID El HU, the crack aperture in the vicinity of the fracture tip can be expressed 
for some S > 0 as: 

w(t, x) = w 0 {t){ 1 - x m ) a ° + u>i(i)(l - x m ) ai + w 2 (t)( 1 - x m ) a2 + o((l - x m ) a2+s ^j , x—*\, (38) 

where powers a,; for respective HF models are given in Table [l] for the case of impermeable rock (qi = 0). 
Note that for the KGD model in toughness dominated regime in the case of perfectly plastic fluid (n = 0), 
the second term of the crack tip asymptotics additionally contains a logarithmic multiplier US¬ 
As a consequence, the particle velocity asymptotics yield: 

v(t,x) = vo(t) + vi(t)(l — x m )P 1 + 0^(1 — x m )P 2 ^J, x—*\. (39) 

The values of /3, are collected in Table [l] 

The value of vq( t) should be greater than zero to make the crack move forward. In the cases under 
consideration vq (t) is defined by no more than the two leading terms of asymptotic expansion for the 
crack opening. 








HF model 

m 

a 0 

Oil 

«2 

Pi 

P2 

Co 

Ci 

C2 

PKN 

i 

1 

n+2 

n+3 

n+2 

2n+5 

n+2 

l 

2 

l 

2 

3 

KGD (viscosity dominated) 

2 

2 

n+2 

n+4 

n+2 

2n+6 

n+2 

l 

2n+2 

n+2 

l 

4 

n+2 

2 

KGD (toughness dominated) 

2 

1 

2 

3—n 

2 

5—2n 

2 

2—n 

2 

l 

2—n 

2 

l 

min(2 - n, ++ 


Table 1: The values of basic constants involved in the asymptotic expansions of w, v and (j> for 0 < n < 1. 


Let us adopt the following symbolic representation: 

lim w n+1 Aw = -C A < 0, (40) 

x —>1 dx L m ~1 ’ v ' 

where C(w) is a functional on the fracture aperture w, related to the form of elasticity operator, while 
Cy i is a known constant. Thus, the formula 


vo 


1 

L m 


1/n 


C A C(w) 


(41) 


is a universal one (valid for both analyzed elasticity operators and all crack propagation regimes) and 
constitutes a relation between the crack propagation speed and the multiplier(s) of the leading ternr(s) 
of the crack opening tip asymptotics. 

The values of Cy i and C{w) for respective HF models are: 

• PKN model 

Ca = ^T2' £M=< +2 , (42) 

• KGD model - viscosity dominated regime 

= < 43 > 

• KGD model - toughness dominated regime 

^ (3 — n)(l — n) mr 

Ca = - 1 -- tan —, C(w)=w 0 +1 wi. (44) 


Taking into account that vq (t) = L'(t), equation (41) can be directly integrated to compute the crack 
length: 


L(t) = 


L 


L +i 


(0) + (— + l) Cjl fc. (' w)dh 


(45) 


Note that, for the toughness dominated regime of the KGD model, the elasticity relation (28) (and 
asymptotic estimate (30)) imply that: 


w 0 (t) = 


(46) 


This, together with relation (44), suggests that in the two limiting cases: small toughness ( I\j —> 0) 
and large toughness (Ki -* oo), the coefficients wq and w i should assume extremely different values, i.e. 
when one of them tends to zero the other magnifies. 
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Combining (41) and (33) with (46), after some algebra and integration, gives the following relation, 
which is valid for any real n: 


L(t) = [L<(0) + tC 


if n V f l I w i + 2 « 

2 U) n ” 


A v 16*7 y 


v (r)tof (r)dr 


(47) 


where £ = K + 2/n+l. For k = 0, the latter coincides with (45) (here m = 2). Taking k = —(n + l)/(2n) 
one can eliminate w o from the formula to have 


m = 


t ¥(„) + l±5 c y»(^ 


W\ n {T)dT 


(48) 


Such a formulation is conducive to computing the so-called small toughness case. However, bearing in 
mind that the value of the functional C(w) is directly evaluated during the computations, relation (45) 
is probably the best choice, and will consequently be used in our algorithm. 

2.3 Reformulation of the problem in terms of reduced particle velocity 

Following |69| . let us introduce a new dependent variable called the reduced particle velocity <fi(t,x): 

x) = v(t, x) — xvo(t). (49) 

The advantages of its utilization have been shown in the recalled paper. The asymptotics of <j> can be 
described qualitatively as (compare Table [I}: 


x) = </>o(i)(l - x m ) c ° + 4>i{t)(l - x m ) 1 ’ 1 + of (1 — x m )‘’ 1+ ' 5 ^, x ->■ 1. 


(50) 


In numerical implementation <f> possesses all the advantages of the particle velocity function. More¬ 
over, an unknown value of the crack propagation speed, Vo, now can be used as an additional regular¬ 
ization parameter in the continuity equation. 

By combining (49) with (22) and substituting the result into the continuity equation (21) one obtains 

1 d 


a modified governing PDE: 


dw 

dt 


Ldx (a4,)+I L W + qi = a ~ 


The boundary condition (36 ft can now be replaced by: 

w(f,O)0(f,O) = q 0 (t), 


while the tip condition (36)2 and the speed equation ( |4l[ ) yield: 

w(t, 1) = 0, (f>(t, 1) = 0. 

Note that the following solvability condition for equation (51) should be satisfied: 

l {jt^ x)+qi ^ x) ) dx+v ml *><!’*)** 


(51) 

(52) 

(53) 

(54) 


and, since w > 0 for any x £ [0,1), it allows one to uniquely define the value of the crack velocity, w 0 , 
at every time step. As a result, equation (51) always has a unique solution with respect to the reduced 
particle velocity </>. 
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Following (22), expressing the pressure derivative in terms of the reduced velocity function gives: 


dp 

dx 


n+1 


W 


{(j) + xv 0 ) r 


(55) 


In this way, the basic system of equations is now built on two dependent variables: the crack opening, 
w(t,x), and the reduced particle velocity, <p(t, x), with an additional relation for the pressure derivative 
I. However, the latter will later be combined with the respective transformed form of the elasticity 


relation (26 1 (271 to eliminate p' x . We will be looking for a solution of the lubrication equation (51), 


under boundary conditions (52), (53) and (35) for the KGD model, initial conditions together with 


the speed equation (33) and resulting expressions for the crack length. Additionally, the fluid balance 


condition (54) and the respective elasticity relation, (24) or (25), together with their inverses, (26) or 


(27), are to be satisfied. 


3 Self-similar formulation of the problem 

Let us assume that a solution to the problem described in the previous section can be expressed in the 
following manner: 

w(t,x) = ip(t)w(x), p(t,x) = L m-u*\ P(x), (56) 


l (ty 


. . . . . '0 l+ ’ l (^)~/ N ,, \ ' ! /’ l+ "(^)l/ \ 

<l(t,x)= L ^ q(x), v(t,x)= L ^ v(x), <j>(t,x)= ^ 4>{x), 

where 7 ) is a smooth continuous function of time. Such a separation of variables enables one 

to reduce the problem to the time-independent form in the case when is a power law or exponential 
function of time. The spatial components of the solution will henceforth be marked by ’hat’-symbol. 
The resulting time-independent formulation is called the self-similar problem. 


3.1 Representation of the solution 

As a consequence of (|56[), the qualitative asymptotic behaviour of the pertinent spatial functions remains 


the same as their time-dependent counterparts (i.e. w(x) complies with (38), v with (39) and <f> with 
(50)). The elasticity equation (23) is converted to: 


p(x) = Aw(x), 

where the self-similar elasticity operator yields respectively: 


(57) 


The inversions of (58) produce: 


w = — - 


TV 


X \ 1 dw r) 

A 2 W(X)= / , 2 2^’ 

Z7r J 0 drj x A — rj z 

(58) 

w(x) =p(x), 

(59) 

K(s , x)ds H— -=Ki\J 1 — x 2 , 
as y 7r 

(60) 


for the PKN and KGD models, respectively. In general, (60) holds provided that the normalized stress 
intensity factor is a function of time defined as: 


K r (t) = kj 


4>(t) 

Vktj 


(61) 
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where: 


K, = — 


p(s) 


Jo Vl — 


ids. 


(62) 


However, as shall be specified later, there is a combination of the self-similar parameters which allows 
one to reduce the problem to the time-independent version for constant Kj. The self-similar particle 
velocity is: 

T dp 1 1/n 


dx 


v{x) = 

with the crack propagation speed computed as: 

vo = (c a jO(w )) " . 


The reduced particle velocity is expressed in the same manner as in (49): 

<p(x) = v(x) — XVq. 

The pressure derivative expressed through <p yields: 

1 


dp 

dx 


„n+l 


XV 0 


)”■ 


Under these assumptions equation (51) can be transformed into the form: 

1 d ~ 

— — (wet)) = -piw - p 2 qi, 

Vo ax 

where p\ and p 2 for two cases of ip(t) are given in Table [2j and: 

qi{t,x) = -if>'(t)qi(x). 

7 

The system has to be equipped with the following boundary conditions 

w(l) = 0, <j>{ 1) = 0, w( 0)^(0) = g 0 , 
with an additional one for the KGD model: 

^1=0. 

dx 1-0 

The fluid balance condition in its self-similar form yields: 

t -< 7 o ~ Pi f wdx - p 2 f qidx = 0. 

Vo Jo Jo 


(63) 


(64) 


(65) 


( 66 ) 


(67) 


( 68 ) 


(69) 


(70) 


(71) 


The foregoing reduction of the problem to the self-similar version is possible for two types of the 
time-dependent functions: 

ipit) = e 7t , ip{t) = (a + f) 7 , (72) 


where a and 7 are some non-negative constants. Values of the respective coefficients in ODE (67), 
together with the pertinent expressions for the crack length for both functions ip, are collected in Table 
[2j Note that the reformulation of the problem to the self-similar form for constant non-zero Kj is only 
feasible for (72)2 with 7 = 71/(2 +n). 
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Type of the self-similar law 

Pi 

P2 

m 

ip(t) = e Jt 

m+n i i 

771+71 

( (m+TpOcA ™+" „ 

2+n 1 1 

( 2 +n )7 

V ( 71 + 2)7 ) 

4>{t) = (a + t) 7 

(m+n) 7 . . 

771+71 

( (m+n)vo \ m +" / , 

n+ 7 ( 2 +rt) 

n+7( n +2) 

[ n+li 2+n)J (V + tj 


Table 2: Values of the auxiliary parameters, pi, p 2 and the crack length, L(t), in the self-similar formu¬ 
lations reflecting different time-dependent behaviour. Here, 0 < n < 1, while mn and Vg are defined in 
Table [I] and (64), respectively. 


3.2 Description of the algorithm for self-similar solution 

The proposed unified algorithm, based on the reduced particle velocity and the speed equation, has a 
common feature for the case n > 0. Namely, equation (64) (and its time-dependent equivalent ©) 
can be directly used to compute the particle velocity. This is no longer true for n = 0, as the equation 
degenerates. Thus the case with n = 0 is a special one and the respective algorithm is built separately, 
although it remains very similar to the basic version. The other limiting case, n = 1, has been completely 
analysed in [B5] using the same approach as adopted in this paper. In fact the only difference, in 
comparison with the general case (0 < n < 1), is that the pressure in the toughness dominated regime 
of the KGD model behaves differently near the crack tip for Newtonian fluid (logarithmic singularity). 
However, as can be easily seen from (44), the speed equation is still valid and has the same form based on 


the coefficients of the asymptotics for the crack opening. Taking all this into account, we built a universal 
algorithm to compute the solution for values of the fluid behaviour index in the range 0.01 < n < 0.99. 
We checked that the difference between the limiting cases (n = 0 and n = 1) and the end of the 
aforementioned n-interval differs negligibly (that is between n = 0.01 and n = 0, as well as n = 0.99 
and n = 1). Nevertheless, while presenting graphical results for the whole domain 0 < n < 1, using our 
numerical and semi-analytical examples, we incorporate the accumulated knowledge about the limiting 
cases as well. All computations were carried out using MATLAB environment on conventional laptop. 

3.2.1 Viscous fluid (0 < n < 1) 

Below we collect a basic set of equations, used later to formulate the universal algorithm for numerically 
simulating hydraulic fractures. The numerical scheme utilizes: 


• the equation defining the reduced velocity following from (671: 

,.i 


i(x) =w)L +p2m ) d( ' 


(73) 


where vg is given in (64), 


inverse elasticity operator (59), (60), written symbolically as: 

l X ) + XVg'j 


• the equation allowing computation of the crack opening, which follows from (66) and the respective 

(74) 


w(x) = B 


the boundary conditions (69), and additionally (70) for the KGD model. 


the solvability condition (71) which follows immediately from (73) and (69)3. 
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For the PKN model operator B in (741 assumes one of the following alternative forms: 



"(v) , 

w n + 1 {rj) V ' 


( 75 ) 


or 


B[v{x)} = 


,1 -I l/(*»+ 2 ) 

(n + 2) / v(rj)dr) 


(76) 


of (28): 


For the KGD variant the operator is directly defined by relation (60) which is a self-similar equivalent 

/i r i 


^K®)] = - 


iv) 


/0 


p,n+l 


0 ?) 


K(j), x)dr) H— y=Ki\/l — x 2 , 


(77) 


where the kernel K(r],x) is defined in (29) and A"/ = 0 for the viscosity dominated regime. 


Note that the general form of the speed equation (33) and the additional boundary condition (35) for 
the KGD model (( [70| in self-similar formulation) are satisfied by the system (73) {77} automatically. 


Remark 2. In our computations for the PKN model we will always use relation (|75j) for the following 
two reasons. First of all, when dealing with the P3D model (see [421 05] ), this is the only form of the 
equation which is available. But more importantly numerical tests revealed that, although formula (76) 


produces more accurate results than (75) in the self-similar formulation, no such advantage exists in 


the transient regime. In this case, the error introduced by the temporal derivative approximation is 


appreciably higher than that resulting from the computation of the operator (74), regardless of its form. 


This feature of the algorithm was also identified in [69]. It is noteworthy that due to this computational 
strategy, the accuracy estimation given in Sectio n |3.4| always describes the upper bound of the error, i.e. 
the accuracy can be improved by using relation ([7617 


3.2.2 Special case: perfectly plastic fluid (n = 0) 


Note that for the perfectly plastic fluid (n = 0) the Poiseulle equation degenerates and formula (63) 
becomes: 

dp 


- w— = 1 , 
dx 


(78) 


which does not allow us to compute the particle velocity. Thus, this case requires special treatment. 

In fact, for the PKN model, one can obtain an analytical solution to the problem provided that the 


leak-off function is predefined and independent of the solution itself. Indeed, combining (78) and ([59]) 
one has: 

w(x) = V 2(1 —x) s, (79) 


which, after substitution into the continuity equation 
particle velocity: 


, allows us to derive a formula for the reduced 


4>{x) = v 0 


1 — x + 


L « da 


2 v/ 27 (l — x) 2 


In the case of the KGD model, relation (60) transforms to a nonlinear integral equation: 


w{x) = - . 

Wo w{i 7) 


[ T7T%^l) + -^A'lt/l - X 2 , 

Jo w(i 7 ) v?r 


(80) 


(81) 


while the reduced particle velocity is computed from (731. 
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Finally, the self-similar crack propagation speed for both the PKN and KGD models has to be 
computed from the balance equation 


vo = 


Qo 


Pi f wdx + p 2 f 0 qidx 


(82) 


3.3 Numerical realization of the algorithm for the self-similar solution 


The solution to the self-similar problem formulated above is sought in the framework of the universal 
algorithm , first introduced in I69j . Even though the algorithm was designed to work for Newtonian 
fluids, in this paper we prove that it is efficient in the case of non-Newtonian fluids as well. 

The adaptation of the numerical scheme addresses new computational challenges related to the 
rheological properties of fluid: i) the order of non-linearity of Poiseulle equation ii) the qualitative 
behaviour of the crack-tip asymptotics dependent on the fluid behaviour index n, iii) the degeneration 
of the Poiseulle equation for perfectly plastic fluid (n = 0). 

The algorithm consists of the following iterative steps: 

• In the first stage, some initial approximation of w = is taker0in a form satisfying the main 

constraints for the crack opening (boundary conditions and asymptotic behaviour). Equations (|7l} 

— (i) — 1 — r 

and (|64|) are utilized to compute v 0 , which is then substituted into (|73|), enabling one to integrate 


the latter and obtain the reduced particle velocity 0W. 


As the governing ODE (671 degenerates at the crack tip, we use the so-called e-regularization 
technique (see e.g. [3D 135] for details), which means that the integration is carried out over the 
truncated spatial interval x € [0,1 — e] , where e is a small parameter. The boundary condition 


(69 >2 is replaced by an approximate one resulting from the asymptotics (50), specified at x = 1 — e. 


The regularized boundary condition is introduced in the form: 

(f>N = Sl(Co, + S 2 (Co, Cl)4>N-2, 


(83) 


where subscripts of </> refer to the indices of nodal points of the spatial mesh (containing N nodes). 
The values of multipliers Si( 2 )(Coi Ci) depend on the particular asymptotic behaviour of function </> 
(see (50) and Table [I]). 

Here £o and Ci are known exponents of the first two asymptotic terms of the function <j >, and as 
such do not change during the iteration process. As a result, functions (f>^ and Vq 1 computed at 
this stage satisfy, together with predefined the: i) fluid balance equation ( pfTj ) , ii) continu¬ 

ity equation (67), iii) regularized boundary condition for ij> ( 83) (equivalent to {69}2), iv) influx 


boundary condition (69)3 - indirectly, through the fluid balance equation. 

At the second stage of each iterative loop, the values of 4>^ and Vq obtained in the previous step 


are utilized to compute the next approximation of using (75) (77). 


While computing the respective integrals in (75) {77}, it is crucial to preserve the appropriate 

asymptotic behaviour of the integrands, resulting from {38} and (50). Hence, computed at 
this stage satisfies the respective elasticity relation and boundary condition (69) i through the 
regularized boundary condition: 


wn — Si (op, a\jWM -i + S2(ao, aq)i/)jv-2- 


(84) 


• The aforementioned two stages of the iterative loop are repeated until all components of the solution 
Vq, 4> and ui have converged to within prescribed tolerances. 


1 Here the superscript j — 1 (j = 1 at the first step) refers to the number of the consequent iteration. 
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Remark 3. Similarly to that found in [55], the performance of the algorithm improves significantly 
when, instead of the dependent variables ui and <fi, one uses the difference between them and their 
leading asymptotic terms: 

Aw = ui — w) 0 (l — x m ) a °, Acj) = <j> — 0o (1 — x m )^°. (85) 


Then the leading terms in the asymptotics of the left and right-hand sides of equations (731 and (75) - 


([77]) are canceled analytically and only the values of Aui and Acj) are iterated. Moreover, while searching 
for the regularization parameter Do, we take into account its relationship with the respective coefficients 
in the solutions asymptotic expansion. This in turn produces a nonlinear equation which is solved using 
the Newton - Raphson method. Finally, the qualitative asymptotic behaviour of the new dependent 
variables Aw and Acj> is known in advance, and the respective exponents should be appropriately adopted 
in the tip asymptotics. Interestingly, sensitivity analysis conducted by the authors proved that, even if 
the employed exponent values are inaccurate, the algorithm remains efficient and stable, however the 
accuracy level deteriorates by one or two orders of magnitude. 

Remark 4. In the case of a perfectly plastic fluid (n = 0) a modified computational strategy has to be 
enforced. For the PKN model it is sufficient to directly employ relations (79) (80), where integration 


of the leak-off function is carried out numerically. For the KGD model, the crack propagation speed Do 
is computed from (82), while the crack opening calculation utilizes (81) with w^~ 1 ' 1 on the right hand 
side. 


3.4 Analysis of the algorithm performance for the self-similar formulation 

Let us now analyze the performance of the algorithm presented above. In the following three subsections 
we investigate the accuracy of computations carried out using the proposed numerical scheme. To this 
end a set of analytical benchmark solutions described in Appendix A is employed. 


3.4.1 Analysis of the algorithm - PKN model 


For the PKN model, we use a benchmark solution in the form (98) for three base functions of the 


type (103). All resulting quantities can be obtained using the methods described in Appendix A. Note 


that this benchmark solution assumes a predefined non-zero leak-off function. All results presented in 
this subsection were computed for 7 = 3 + 2n > whidb in the applied scaling, corresponds to a constant 
injection flux rate. 

In this subsection we use various error measures in order to demonstrate the algorithms peculiarities 
and to define an optimal set of accuracy parameters, used later on to investigate the performance of 
the proposed numerical scheme. We will consider both the average and maximal relative errors of the 
following: the crack opening, Sw, the particle velocity, 5v, and the reduced particle velocity, Scf>. The 
relative errors over x are defined by the following standard form: 


5m av 



VJ]Sl) 2 dx ■ 



5m 


max 


max 

0<x<l 


\m(x) — wn{x) I 


\m{x)\ 


( 86 ) 


where m stands for the benchmark solution, while wn refers to the respective numerical result. 

In Fig. [T] the average error of the basic variables, as functions of the mesh density, N, and the 
fluid behaviour index, n, are depicted. The errors were estimated for a number of nodal points varying 
from 10 to 200 , where the spatial mesh density was increased at both ends of the interval, in a manner 
described in [55]. The smallest computed value of n was moved away from zero (to 0.01), as for the 
PKN model the solution for the perfectly plastic fluid is found in the framework of a simplified scheme 
- see Remark[4](it can, in fact, be derived analytically for the considered benchmark). 
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Figure 1: PKN model - average relative error in the self-similar solution for different numbers of nodal 
points, IV, and values of the fluid behaviour index, n: a) the crack opening w, b) the particle velocity v, 
c) the reduced particle velocity 



Figure 2: PKN model - maximal relative error in the self-similar solution for different numbers of nodal 
points, N, and values of the fluid behaviour index, n : a) the crack opening w, b) the particle velocity v, 
c) the reduced particle velocity </>. 


From Fig. [l] it follows that the average computational error is hardly sensitive to the value of fluid 
behaviour index, n. Indeed, it is only for the largest number of nodal points that a discernible difference 
in accuracy for different n can be observed. Simultaneously, for any fixed n, a constant improvement in 
the solution accuracy when increasing the number of nodal points N can be observed over almost the 
whole range of N. Some irregularities (local minima) in the error distribution can be observed for the 
particle velocity and the reduced particle velocity for N around 150. For a fluid behaviour index greater 
than approximately 0.5 when taking N > 150 the error is no longer reduced. No noticeable difference 
between the error distributions for v and <fi is observed. 

In Fig. [ 2 ] we depict the distributions of the maximal relative errors for w, v and 4> as functions of N 
and n. It shows that the general trends in accuracy (low sensitivity to the value of ?z, error reduction 
with growing N) are similar to those reported above. However, the level of the respective errors for the 
crack opening and the reduced particle velocity are up to two orders of magnitude higher than they were 
for the average error. This trend does not hold for the particle velocity, where the maximal error is of a 
similar order to that for the average formulation. 

In order to explain the difference between the average and maximal errors let us analyze the error 
distribution over x for a given number of nodal points, N = 100. In Fig. [d]we present the relative errors 
of w, <j> and v as functions of x and n. The error distribution of the crack aperture and the reduced 
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particle velocity is rather non-uniform, with a distinct increase in the near-tip region. In the case of 
the particle velocity (Fig. [3 }d) one does not observe the accuracy deterioration at the crack tip, which 
is a result of a different asymptotic behaviour. Indeed, the spike growth of Sw(x) and S(f>(x) is a direct 
consequence of the fact that both variables vanish at the crack tip (compare with (69)i and ( 69 ) 2 ) while 
the particle velocity tends to a non-zero value 0(1) = Vo 7 ^ 0. Moreover, the graph of Sv indicates the 
quality of computations for the crack propagation speed, Vo, itself. The reduced particle velocity, (f>, and 
particle velocity, 0, exhibit some tendency towards error reduction at the crack inlet for lower values of 
n, however this does not result in a decrease of the maximal error, which is similar for the whole range of 
n. Thus, the relative error magnification at the crack tip has little influence on the average error values, 
as the accuracy over most of the x interval is up to three orders of magnitude better. 



Figure 3: PKN model - distribution of the relative error for the self-similar solution: a) the crack opening 
w, b) the particle velocity v, c) the reduced particle velocity <j). The mesh is composed of TV = 100 nodal 
points. 


As a consequence of the foregoing analysis, from now on we will present as graphical information only 
the following measures of the solution accuracy: i) the average relative error over x of the crack opening 
and the particle velocity; ii) the relative error of the crack opening and the particle velocity for a fixed 
value of TV. The first group of parameters is very useful when analyzing the general trends in accuracy 
for various TV and n. Moreover, it has been shown that the average error for v is practically the same 
as that for <j>, which in light of the straightforward physical interpretation of the former makes it more 
convenient for use in accuracy analysis. The second group of parameters (point-wise relative errors) 
provide an illustration of the error distributions local peculiarities, while simultaneously depicting the 
solutions maximal errors. 

Before we move on to the KGD models, we would like to note that the values of computational 
parameters such as the regularization parameter e = 10~ 6 and mesh density parameters defining the 
distribution of points on the spatial interval (0,1) were kept identical in all computations. It is clear 
from the results presented above that, by optimizing these parameters for each specific value of the fluid 
behaviour index, n, and number of nodal points, TV, one can reduce the error further, and there is an 
apparent potential for greater accuracy improvement with TV growing beyond 200. However, the level of 
computational accuracy reached already is at least a few orders of magnitude better than any reported 


in the literature (see results in Section 3.5). 


3.4.2 Analysis of the algorithm - KGD model in viscosity dominated regime 

Let us now perform accuracy analysis for the viscosity dominated regime of the KGD model. The 


benchmark solution is based on representation (98) for three test functions of type (105). All resulting 


quantities can be obtained in the manner described in Appendix A. As for the PKN model, this bench- 
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mark solution assumes a predefined non-zero leak-off function. All results presented in this section were 
computed for 7 = ^77, which corresponds to a constant influx rate. 

We begin by investigating the relationship between the accuracy of computations, mesh density 
(number of nodal points N) and the fluid behaviour index, n. For different numbers of nodal points N 
(ranging from 30 to 200) the computations were performed for n continuously changing between 0 and 1. 
Again, the mesh density was increased at both ends of the spatial interval. The lower limit of N was set 
to 30 instead of 10, as proper numerical computation of the inverse elasticity operator (60) necessitates 
finer meshing than its equivalent in the PKN model (identity operator). It should be emphasized that 
it is still possible to obtain very accurate results even for N < 30, however this requires modifying the 
nodal spacing at the ends of the interval for each value of n. 

Similarly to that shown in Fig. [T] for the PKN model, the numerical simulations performance in 
the KGD viscosity dominated case is demonstrated using the average relative error over x for the crack 
opening, dw av > and the particle velocity, 5v av , with the results presented in Fig. [IJ It shows that the 
level of both errors is the same, however they are distinctly greater (up to two orders of magnitude for 
large N) than was the case for the PKN model. In the analyzed range of N, increasing the mesh density 
leads to a monotonous reduction in the error of computations. However, for N < 200 the stabilization 
level is still not achieved. Low sensitivity of the results to the value of n is observed. Only for coarse 
meshing does the solution exhibits some susceptibility to this parameter. The trends reported for 5(j) av 
(not depicted here) were the same as those for Sv av . 



Figure 4: KGD model (viscosity dominated regime) - average relative error for the self-similar solution 
with different numbers of nodal points, N, and values of the fluid behaviour index, n: a) the crack 
opening ut, b) the particle velocity v. 

As previously for the PKN model, here we also present the spatial distribution of the solution error 
for N = 100. The analyzed measures of error are: the relative error of the crack opening, 5w, and the 
relative error of the particle velocity, 5v. Additionally, the absolute error of the net fluid pressure, A/3, is 
introduced, as later on we are going to compare our results for the pressure with others that are available 
in literature. Note that in the case of the PKN model, the error of the net fluid pressure was solely 
defined by the error of the crack opening. We do not show the pressures relative error, as the function 
p assumes a zero value at some point in the interval. All computations were carried out for N = 100 
nodal points and non uniform spatial meshing (the same for all n). The respective results are presented 
in Fig. [5| 

First and foremost, regardless of the error measure, the same accuracy level is retained over the 
whole interval of n. This proves that the proposed numerical scheme is highly stable. The relative error 
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Figure 5: KGD model (viscosity dominated regime) - error distribution for the self-similar solution: a) 
relative error of the crack opening 5ui , b) relative error of the particle velocity dv, c) absolute error of 
the net fluid pressure A p. Mesh is composed of N = 100 nodes. 


of the crack opening does not exhibit a magnification at the crack tip, which was evident in the PKN 
model. The maximal value of Sw remains similar to that in the PKN case (compare Fig. [3^i), but its 
location has moved away from the crack tip. In general, the error level over most of the spatial interval 
(except for the near-tip region) has increased in comparison with the PKN model by at least an order 
of magnitude for the crack width and the particle velocity. This is not a surprise as the KGD problem 
is more complicated than the PKN model. Surprisingly, the accuracy of the net fluid pressure is of the 
same level even near the fracture front, where p exhibits singular behaviour. 


3.4.3 Analysis of the algorithm - KGD model in toughness dominated regime 


The benchmark solution, utilized here for accuracy analysis, is based on representation (981 for five base 


functions of type (|105j and ( |110[ ). The applied power 7 = corresponds to a constant injection flux 


rate and simultaneously allows the reduction of the time-dependent problem to its self-similar version 
for a constant (independent of time) value of the normalized stress intensity factor. The self-similar SIF 
for this benchmark was Kj ss 1 and depended on n. 

Let us start by analyzing the average relative errors, similar to that conducted previously for the 
PKN model and the KGD model in viscosity dominated regime. In Fig. [6] we depict the average error 
of the crack opening and the particle velocity for various N (30 < N < 200), with n in the range of 
interest. As previously, the adaptive spatial mesh (the same refined meshing at the ends of the interval 
for all values of n) was used. 

The results show that again the accuracy of computations is not sensitive to the value of the fluid 
behaviour index n. The error distributions for both analyzed dependent variables are very similar. The 
error level is slightly higher than it was for the viscosity dominated regime of the KGD model (compare 
Fig. [4|. However, it is sufficient to take as few as 50 nodal points to have the average relative error 
of order 10^ 5 . In the analyzed range of N an increase in the mesh density produces a constant error 
reduction. No saturation level is achieved for N < 200. 

In the next stage of our analysis we investigate the spatial distribution of relative errors. In Fig. [7] 
we present the relative error for w and v and the absolute error of p as functions of x and n. In general, 
the accuracy for the crack opening and the particle velocity is about half an order worse than it was 
for the viscosity dominated regime and does not change with n. No error magnification at the crack tip 
is observed. Also, the error of the net fluid pressure proves to be very stable for various n. A small 
increase in A p can be seen at the crack tip for growing n, however considering the tip singularity of p 
this phenomenon does not detract from the high quality of the results. 
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Figure 6: KGD model (toughness dominated regime) - average relative error for the self-similar solution 
for different numbers of nodal points, N, and values of the fluid behaviour index, n: a) the crack opening 
w, b) the particle velocity v. 



Figure 7: KGD model (toughness driven regime) - error distribution for the self-similar solution: a) 
relative error of the crack opening 5w 1 b) relative error of the particle velocity 6v, c) absolute error of 
the net fluid pressure A p. The mesh is composed of N = 100 nodes. 


3.4.4 Summary of the presented results 

In order to present a brief summary of the foregoing accuracy analysis let us collect representative values 
of the computational error for the main components of the self-similar solution: the crack opening, ui, 
and the particle velocity, v. The results are given in Table [3] For each hydraulic fracture model and 
accuracy parameter, the presented figures should be considered as averaged over n. In some cases the 
difference between extreme values for different n was negligible, in others there was a distinct discrepancy 
but they still remained within the same order of accuracy. The quality of computations for the fluid 
flow rate can be easily reconstructed from the presented data (bearing in mind that q = vjv). One can 
also estimate the accuracy of the net fluid pressure for each HF model under consideration from Fig. 

Fig. §: and Fig. [7}:. All computations were carried out using a mesh composed of N = 100 nodal 
points, with increased density at the ends of the interval. 

To enrich the presented investigation even further, we have also conducted additional tests for differ¬ 
ent values of the self-similar SIF, Kj = 100,1,0.1,0.01. To this end, the analytical benchmark solutions 
from Appendix A have been used for different magnitudes of parameters A j. We wish to underline that, 
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even for the lower values of AT/, the problem should not necessarily be viewed as the so-called ’small 
toughness’ regime. It results from the presence of a non-zero leak-off, which in the considered class of 
benchmark solutions is properly combined with the fracture opening and the particle velocity to balance 
the continuity equation. Thus it can alleviate, to some degree, the computational difficulties created by 
this extreme regime of crack propagation. Indeed, the results presented in Table [3] seem to confirm this 
prediction, as only a limited deterioration in the accuracy is observed for the two lowest values of I\j. 
The problem of small toughness will be analyzed in the next subsection. 


Variable 

PKN 

KGD 

KGD toughness dominated 

viscosity dominated 

Kj = 100 

K, = 1 

kj = o.i 

kj = o.oi 

6w av 

10 -9 

3•10~ 7 

5•10~ 7 

10" 6 

6 •10 -6 

10' 5 


10 -9 

io - 7 

IO" 7 

5•10 -7 

10" 6 

10~ 6 

SvJmax 

4-IO" 7 

4- 10" 7 

5•IO’ 7 

3•10" 6 

10" 6 

4-10" 6 

fi'Vmax 

10" 8 

2 •10~ 7 

10" 6 

5•10 -6 

IO" 5 

10~ 5 


Table 3: The level of computational error for the self-similar crack opening, w and the particle velocity, 
v, computed by the universal algorithm for different HF models, based on analytical benchmarks from 
Appendix A. For all computations the non-uniform mesh, distributed over the integral (0,1), consisted 
of IV = 100 points and was more dense near the endpoints. Meshing parameters, and all other technical 
parameters, remained identical irrespective of the chosen model. 


3.4.5 KGD model: small toughness case 

It is well known that the small toughness case constitutes an extremely challenging computational 
problem due to the immense reduction of the process zone [55] . It necessitates special measures to 
stabilize the computations and preserve the solution accuracy, which results in an appreciable increase 
in the computational cost (in comparison with a moderate toughness case, the solution time was 100 
times higher). In our test we shall utilize a reference solution for the viscosity dominated regime of the 
KGD model, namely a numerical solution obtained for qi = 0 and go = 1. It is obvious that if for zero 
leak-off one imposes the same influx value in the toughness dominated regime, then for sufficiently small 
Kj the solution should approach that of the viscosity dominated regime. 

The strategy for our test is as follows. For three values of the fluid behaviour index n = {0, 0.5,0.9} 
we conduct computations while steadily decreasing the self-similar stress intensity factor, A'/. At every 
step we investigate the relative deviation of the crack opening in the toughness dominated regime (wt) 
from the viscosity dominated reference solution (wy): Sw = \(wt ~ Wy)/wy\. The results for Kj = 
{0.1,0.05,0.01} are depicted in Fig. [8j 

It shows that for a fixed AT/ the level of relative deviation Sw depends on the fluid behaviour index, 
n. In general the lower the value of n being analyzed, the bigger the reduction in Kj that is needed to 
achieve a certain level of discrepancy. Naturally, the biggest deviation between wt and wy is observed 
at the crack tip, where the respective asymptotics do not match. However, it is sufficient to take 
Kj = 0.05 to obtain a relative difference from the viscosity dominated solution of the order 10” 3 for 
any n on practically the whole crack length. One can conclude that, when Kj < 0.05, it is acceptable 
to use the solution of the viscosity dominated variant of the problem in order to achieve the accuracy 
needed in any practical application. These results provide an independent indirect verification of the 
computational accuracy of the algorithm. However, accurate computing in the limiting case of small 
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Figure 8: Crack opening: the relative deviation of the numerical solution for the small toughness KGD 
model from the KGD viscosity dominated reference solution for different values of the self-similar stress 
intensity factor, Kj. The mesh is composed of N = 100 nodes. For each I\j appropriate value of the 
regularization parameter e was taken. 


toughness necessitates proper handling of the computational parameters (the regularization parameter 
e and mesh density). 


3.5 Comparison with other available results 

Having identified the level of accuracy achievable by the presented numerical scheme we can now conduct 
computations for the reference data available in literature. In the following sections we shall present a 
comparative analysis for a few solutions delivered by other authors for PKN and KGD models. 


3.5.1 PKN model 


In |41| an analytical solution for the PKN model was delivered in terms of an infinite power series for 
the so-called ’’proper variable” w n+2 and the particle velocity v. Due to the complexity of the problem 
the coefficients of the series are represented in the form of recurrent sums, which makes computing them 
a complicated task in and of itself. For this reason, the author presented explicit formulae for the first 
five coefficients in respective expansions. Moreover, the representation of the solution also contains an 
additional parameter, the self-similar crack length (see [41] h £*, which is to be determined from the 
influx boundary condition. Unfortunately, the slowest power series convergence over the whole range 
of x takes place for x = 0. Thus one cannot expect good accuracy for the self-similar crack length 
approximation when using merely five terms. On the other hand, the author observed that the unknown 
£* lies between £* = 1.04004 (for n = 0) and £* = 1.00101 (for n = 1). This led him to conclude 
that the averaged value of the self-similar crack length £* = 1.02 is sufficient to produce, with 5-term 
approximation, a solution whose error does not exceed a few percent for any n. 

In the following, with the accurate computational results obtained using the universal algorithm, we 
shall: i) recreate the respective numerical values of the self-similar crack length £* pertaining to the 
problem formulated in m - here an approximate formula for £*(n) will be given; ii) estimate the error 
of the 5-term solution approximation, given in mi for three variants of £*: the actual (numerical) values, 
the averaged value = 1.02, and a linear approximation between the extreme values, £*(0) = 1.04004 
and £*(1) = 1.00101, given in [21] . 

All computations were carried out using a mesh composed of 200 points which, according to the 
accuracy analysis given in Section 3.4.1 provides a solution accuracy of up to order 10~ 9 for the crack 
opening and 10~ 8 for the particle velocity. Appropriate rescaling of the problem to the formulation used 
in [44] was carried out under the assumption, accepted there, that the self-similar crack propagation 
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speed is a known function of £*. A number of computations for different values of the fluid behaviour 
index n allowed us to recreate the relation £*(n) numerically. Based on the numerical simulations we 
have established an approximate formula for the self-similar crack length £*(n): 


2n+3 



(1 — n ) 


2-\/3 0.07221n - 0.04222 

- 4 - Tl - 

3 n 2 + 3.247 n + 2.458 


+ 0.99831387n, 


(87) 


which imitates the accurate numerical results with an accuracy of order 10“ 7 . 

Having the values of £*(ro) for an arbitrary fluid behaviour index enables us to assess the accuracy 
of the 5-term solution approximation proposed in |41j . We have checked that the error of the 5-term 
approximate solution remains exactly the same when recreated using relation (87) or with the accurate 
numerical value of £*. The respective results are presented in Fig. [9] 

It shows that for both dependent variables such an approximation provides a maximal error for w 
of the level 10 -5 and 10~ 4 for v. The accuracy improves with decreasing n. The crack opening, w, has 
greater potential for further error reduction. 


a ) 


b) 


5w 



i 


Sv 



i 


Figure 9: PKN model comparison of the 5-ternr approximation solution presented in [ IT] with the 
accurate numerical solution for N = 200: relative errors of: a) the crack opening, w(x), b) the particle 


velocity, v(x). Self-similar crack length was computed according to (87). 


For comparison, in Fig. [l0]we depict the accuracy of the 5-term approximations for w and v, obtained 
using the averaged value of the self-similar crack opening £* = 1.02 suggested in m- Now, Sw and 5i> 
are of a similar level with maximal values (around 1 percent as predicted) located at n = 0 and n = 1. 
The best results are obtained in the vicinity of n = 1/2, where the averaged £* is close to the real value. 
A natural suggestion emerging from this test is that a linear approximation of the self-similar crack 
length, based on the two extreme values (for n = 0 and n = 1), should produce better results than the 
constant averaged value. In Fig. EH we depict the errors of w and v obtained for such a variant. The 
maximal errors have been reduced by one order of magnitude with respect to the previous case. The 
best coincidence with our numerical solution is obtained at the ends of the interval n. 

In order to cross-check the credibility of the respective results, let us investigate to what degree they 
satisfy the influx boundary condition. In the problem considered in [4lj a predefined unit value of the 
influx was assumed (the normalization of the governing equations taken as in [411 ~): qo = 1. In Fig. 12 
we show the relative errors of <70, recreated as a product of the crack opening and the particle velocity 
(qo = w(0)v(0)), for our numerical solution and different variants of the 5-term approximation. As can 
be seen in the figure, the numerical solution fulfills the influx boundary condition with a stable accuracy 
of the level 10 -8 over the whole range of n. The 5-term approximation with £* computed from (87) yields 
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Figure 10: PKN model - comparison of the 5-term approximation solution presented in m with the 
accurate numerical solution for N = 200: relative errors of: a) the crack opening, w(x), b) the particle 
velocity, v(x). The averaged (constant over n) self-similar crack length was taken: £* = 1.02 (see page 
16 in [41|). 


a ) 


b) 


Sw 



i 


6v 



i 


Figure 11: PKN model - comparison of the 5-term approximation solution presented in m with the 
accurate numerical solution for N = 200: relative errors of: a) the crack opening, w(: r), b) the particle 
velocity, v(x). The self-similar crack length, £*, was taken as a linear function of n. 


more accurate results for lower values of the fluid behaviour index, with a monotonous error increase 
up to a level of ICR 4 for growing n. As anticipated, for £* linearly dependent on n one obtains the best 
accuracy for the extreme values of the latter, while the averaged £* = 1.02 produces optimal results for 
n « 0.5. To summarize, the 5-term approximation given in m, complemented with the approximated 
value of the self-similar crack opening £*(n) provided by formula (87), can be considered a very good 
benchmark for the PKN problem with non-Newtonian fluid. 

Recently, two new analytical solutions of the PKN problem have been proposed in m for imperme¬ 
able rock. One of them refers to the case when the fracture volume does not change with time for a zero 
injection flux rate (shut-in regime). We will not analyze this benchmark as it has no relevance to the 
real (practical) situation. However, it demonstrates that there exists a set of nontrivial solutions to the 
problem if the initial crack opening is different from zero. 

The second benchmark solution proposed in m assumes that the particle velocity is constant along 
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Figure 12: PKN model - fulfillment of the influx boundary condition by: cq - the numerical solution, 82 
- the 5-term approximation with £* computed from (87), ^3 - the 5-term approximation with averaged 
£* = 1.02, £4 - the 5-term approximation with £* linearly dependent on n. 


the entire fracture: v(t,x) = Vq( t) = L'(t). Below we present a reproduction of this solution in terms of 
the self-similar formulation of type (72)2 for 7 = l/(n + 2). It can be proved that the following set of 
equations constitutes a solution to the considered self-similar problem: 


w(x) = t&o(l — :r) 1/, (" +2 \ wo = (n + 2) ” (n + 2 > 


( 88 ) 

(89) 


v = v 0 = (n + 2) . 

Note that this solution can be directly obtained from the analytical solution given in m- 

Having this benchmark we can now conduct a similar accuracy analysis to that presented in the 
previous section. In Fig. [13] - [14] we depict the average relative error, over x , of the crack opening and 
the particle velocity for 10 < N < 200. It shows that the error level is hardly sensitive to the value 
of n. In the analyzed range of N no accuracy saturation is observed. Although the general trends in 
error distribution are similar to those reported for the non-zero leak-off benchmark in Fig. [T] the overall 
error level is slightly higher. This allows one to conclude that the accuracy of the computations by the 
universal algorithm is related to the type of problem being analyzed. We believe that this feature is not 
unique to the proposed scheme, as similar trends were reported in [3% SB I, :55]. 


Fig. 14 illustrates the distribution of the relative error of w and v over the spatial interval for N = 100 
with different values of the fluid behaviour index n. For all analyzed n we see an almost constant error 
over the whole range of x. No error magnification at the crack tip (as in Fig. [3]) is observed. Indeed, 
in this case Swq and 8vq define the error level over x (compare (88) and (89)). Compared to the results 


presented in Fig. [3| the level of 8v remains the same while 8ui is reduced by one order of magnitude. 
Interestingly, a pronounced reduction in 8w takes place for n approaching zero. We believe that this 
fact further confirms the good performance and stability of our scheme, as the analyzed benchmark 
degenerates for the perfectly plastic fluid. 


3.5.2 KGD model in viscosity dominated regime 

Let us now consider a reference solution given in [Tj. The authors analyzed the viscosity dominated 
regime of the KGD model for a number of shear-thinning fluids. They also assumed a constant influx 
and impermeability of the rock formation (qi = 0). Semi-analytical solutions to the problem were 
proposed for a variety of shear-thinning fluids (n = 0, 0.1,0.3, 0.5,0.7, 0.9,1). 
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Figure 13: PKN model - verification of accuracy of the computations by the universal algorithm against 
the analytical benchmark (881 (89): average relative error of the numerical solution with different 
numbers of nodal points, N, and values of the fluid behaviour index, n: a) the crack opening w, b) the 
particle velocity v. 



Figure 14: PKN model - verification of accuracy of the computations by the universal algorithm against 
the analytical benchmark (88) - (89) for TV = 100: average relative error of: a) the crack opening w, b) 
the particle velocity v. 


In the following we compare our numerical results with those given in the recalled paper in terms of: 
i) self-similar crack opening, ii) self-similar fluid pressure, iii) self-similar fluid flow rate, iv) self-similar 
particle velocity. Although in [TJ there is no data for the particle velocity, it can be easily obtained 
through the fluid flow rate and the fracture opening (v = q/w ). As for results obtained using the 
universal algorithm, the particle velocity is retrieved directly from the reduced particle velocity, while 
computing the fluid pressure requires additional postprocessing (integration). The predefined influx 
value was taken from the data given in [T] (p.591 Table I). 

Our computations were carried out on an adaptive mesh composed of 200 nodes, which according 
to the conducted accuracy analysis should produce errors of order 10 -7 — 10 -8 . Comparison of the 


numerical results with those from |l] is shown in Fig. 15 - 17 
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As there are two ways of calculating the fluid flow rate q in [T] - either by a sum of the base functions 
(U/ m o in equation (41) in the aforementioned paper) or using q = wv (where v is computed from (63)), we 
shall analyze both. From now on we will refer to these methods as q^ and (f 2 \ respectively. Pertaining 
formulae for the particle velocity are defined as: = q^/w, while v^ is computed according to (63). 

Note that for the exact solution obtained using both computational methods should produce the same 
result. For the perfectly plastic fluid (n = 0) neither q^ nor v^ can be obtained, as equation (63) 
degenerates to yield (78) and thus it can no longer be used to define the fluid flow rate q or the particle 
velocity v as it was explained in Section 3.2.2[ Therefore the respective curves are not presented in Fig. 
[Mr and [Mr. 




Figure 15: KGD model (viscosity dominated regime) - comparison of the semi-analytical solution pre¬ 
sented in [T] with our numerical solution for N = 100: a) the crack opening w, b) the fluid pressure 

P- 


As can be seen for the crack opening and fluid pressure, the curves corresponding to the universal 
algorithm and to the solution obtained by [1] are indistinguishable from one another in the presented 
scale (the x interval in Fig. 15 r was truncated because the pressure tends to infinity as x —> 1). When 
one analyzes the particle velocity it turns out that the values given by the universal algorithm are 
in good agreement with v^ over almost the entire interval. A distinct divergence in results can be 
observed in the near-tip region, where the quality of approximation v^ apparently deteriorates. On the 
other hand, when retrieving the particle velocity from the data given in [T] as v^\ we have very good 
coincidence in the vicinity of the fracture front, with the results diverging from each other as the crack 
inlet is approached. In general, better coincidence of with our accurate numerical result for v takes 
place for lower values of fluid behaviour index n, while a reverse relation is observed for v^ 2 \ A similar 
situation to that reported for v has been revealed for the fluid flow rate. Again, q^ gives better results 
away from the crack tip, while q^ produces better coincidence in the near-tip region. 

The above analysis confirms the credibility of our solution, which together with the previous accuracy 
estimation allows us to now treat it as the numerical benchmark data. Following the idea from [Ti and 
(40] . we propose a new improved approximation of the dependent variables analyzed above. It provides 
a solution to the problem with far greater accuracy than any other known semi-analytical formulae and 
thus can be treated as the reference data itself when testing computational algorithms for hydraulic 
fractures. The new solution approximation presented below is valid for any value of the fluid behaviour 
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Figure 16: KGD model (viscosity dominated regime) - comparison of the semi-analytical solution pre¬ 
sented in [T] with our numerical solution for N = 100: the fluid flow rate q. 




Figure 17: KGD model (viscosity dominated regime) - comparison of the semi-analytical solution pre¬ 
sented in [lj with our numerical solution for N = 100: the particle velocity v. 


index. Let the fracture opening, w, the fluid pressure, p, and the particle velocity, v, be expressed as: 


w 


( x ) = w 0 (l - x 2 )"+ 2 + W\ (1 - z 2 ) 4/3 + 


w 2 


\/l — x 2 — \ (1 — x 2 ) 3 / 2 — x 2 In ^ ^ - 


(90) 


p{x) 


a) /i 

(n + 2)7r 2 ^' 1 \2 


2 

n + 2 ’ 


1 ; 1 / 2 ; x 
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+ ~ X ) 
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v(x) 
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2=1 


(91) 

(92) 
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while the respective approximation of the fluid flow rate, q , can be easily obtained from the product of 


the fracture opening (90) and the particle velocity (921. 


In the equations above, B is the beta function, 2 -F 1 denotes the Gauss hypergeometric function, C _4 


is as expressed in (431 and wq is computed using the formula (0 < n < 1): 


wo = 


n + 2 


1 \ n/(n+ 2) 

t \ l/(n+2) 


c: 


(93) 


Coefficients Pi(n) and Vi(n) are functions of the fluid behaviour index n. Their values are collected 

in Table |4] (Appendix B). 


In Fig. 18 19 we present a comparison of the new approximation (901 (921 and the semi-analytical 


formulae from fl] with the exact numerical solution. Bearing in mind the accuracy level provided by the 
universal algorithm we claim that the depicted deviations describe, in fact, the error of the respective 
approximations. Fig. |18f r illustrates the relative deviations of the crack opening computed as the 
averaged, S av , and maximal, S max , values. For the net fluid pressure (Fig. 18 :>) we show only the 


absolute deviations (again in the mean and maximal formulations), as the pressure graph intersects the 
x axis. The deviation of the fluid flow rate and particle velocity are presented in Fig. [19] Here, the data 
from [I] was computed in two alternative ways: as gW, and q( 2 \ i)( 2 \ 
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Figure 18: KGD model (viscosity dominated regime) - comparison of the new approximate solution (901 


(92) and the semi-analytical solution from [T| with accurate numerical results for N = 100: a) the 
average and maximal relative error of the formulas for the crack opening 5w av and Sw max ; b) average 
and maximal error of the fluid pressure A p av and A p max . 

The first straightforward conclusion that emerges from the presented results is that, regardless of the 
error measure and analyzed component of the solution, the new approximation gives at least one order 
of magnitude higher accuracy than the approximation from pQ. Moreover, better stability of the results 
produced by (90) - (92) is observed when varying n. When analyzing the graphs for Sq and Sv it shows 


that the quality of the solution from [I] depends essentially on the way in which the fluid flow rate and 
the particle velocity are computed. In general, the method referred to as q( 2 \ v^ is better over almost 
the entire range of n, while both approaches are of similar accuracy for n close to unity. However, this 
conclusion does not reflect the total complexity of the problem. The spatial distribution of respective 
errors reveals that q^ and v^ are much better than q ^ and over almost the whole crack length, 
except for the near-tip region where pronounced error magnification takes place. On the other hand, 
g*' 1 ) and prove to be very accurate when approaching the crack tip and diverge from accurate results 
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Figure 19: KGD model (viscosity dominated regime) - comparison of the new approximate solution (901 
(92) and the semi-analytical solution from [T] with accurate numerical results for N = 100: a) the fluid 
flow rate q, b) particle velocity v. Superscript here defines which evaluation formula for the results from 
PQ has been employed. Both average relative, 5f av , and maximal relative, Sf max , errors are presented. 


with decreasing x. These trends can be observed in Fig. [16]- Fig. [T7] Moreover, exactly the same trend 
was identified for Newtonian fluid in [BUj . However, if one needs to utilize the approximation given in [T] 
it seems that the optimal way to do so should be based on appropriate merging of both representations 
for q and v. 

In order to supplement the foregoing analysis we shall now present a detailed description of the 
accuracy of approximation (90) (92). In Fig. 20 - Fig. 21 we depict the error distributions over n and 


x for the crack opening, Sw, the net fluid pressure, A p, the fluid flow rate, Sq, and the particle velocity, 
5v. As previously, the relative errors are shown for ui, q and v, while inaccuracy of p is described by its 
absolute measure. The first conclusion from the presented data is that the new improved approximation 
provides accurate and stable results regardless of the value of the fluid behaviour index, n. None of the 
presented errors exhibit sharp magnification at the crack tip, which means that the proposed formulae 
preserve the tip asymptotics very well. The maximal relative errors of the crack opening and the fluid 
flow rate are far below one percent (and far below the errors generated by the solution from [Tj), while the 
error of the particle velocity is even lower and does not exceed 2 • 10 -4 . The net fluid pressure gives very 
accurate results even at the crack tip, where the singular behaviour holds. Furthermore, discrepancies 
between the improved approximation and the reference numerical data decrease as the fracture front is 
approached. 


In the light of the foregoing analysis we can conclude that the new improved approximation (90) 


- (92) constitutes a very good imitation of the reference solution for any 0 < n < 1 and thus can be 
itself used as a semi-analytical benchmark for testing other numerical algorithms. In any case, it is 
appreciably more accurate than the solution proposed in .1] and, unlike the former, avoids rapid error 
growth at the crack tip. 


Remark 5. Respective formulae from (90) (92), as defined on the basis of a pre-computed numerical 


solution, should be considered independent from each other. This means that if one wants to recreate 
e.g. the particle velocity by substituting (90) and (91) to (63), the resulting accuracy may be lower than 


that reported for (92). 
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Figure 20: KGD model (viscosity dominated regime) - comparison of the new approximate solution (901 
(92) with accurate numerical results for N = 100: a) relative deviation of crack width ui, b) absolute 
deviation of fluid pressure p. 



b) 


6v 



Figure 21: KGD model (viscosity dominated regime) - comparison of the new approximate solution (90) 
- (92) with accurate numerical results for N = 100: relative deviation of: a) fluid flow rate q , b) particle 
velocity v. 
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3.5.3 KGD model in toughness dominated regime 


Let us now compare the numerical solution produced by the universal algorithm with the reference data 
for the toughness dominated regime given in m There the author presents a semi-analytical solution 
for the crack opening and the net fluid pressure in terms of a series expansion for three values of the 
fluid behaviour index: n = {0,0.5,1}. Results for n = 1 have already been analyzed in {B5J. Now we 
shall investigate the remaining variants (n = 0 and n = 0.5). Following the conclusions from Section 


3.4.3 we conducted the computations using a spatial mesh consisting of N = 100 nodes, which previous 
accuracy analysis showed provides a solution error of at most order 10 _o . 

The plots of the crack opening, vj, and the net fluid pressure, p, for n = 0 and n = 0.5, together with 
respective data from HU, are drawn in Fig. |22[ Since the normalized values of the fracture apertures 


for these cases are very similar (even near the crack tip), they are presented in the separate graphs 
(Fig. 22 1 and Fig. j22jr). Interestingly, the pertaining curves for the net fluid pressure are not that 
close to each other, especially in the near tip region. An explanation for this phenomenon is as follows. 
Although the leading terms of w are of the same order (a 0 = 1/2, see Table [l]), the second terms 
depend on the fluid behavior index (aq = (3 — n)/2). Consequently, on substitution of the asymptotic 


representation of w into the elasticity operator (58)2, the square root term produces a constant value, 


while the second term gives an expression of power law singularity that depends on the fluid behaviour 
index (p(x) ~ (1 — a, 2 )! 1 -™)/ 2 as x —> 1). For this reason the tip behaviour of the net fluid pressure is 
essentially different for both cases. Moreover, it hinges on the value of the multiplier of the second term 
of asymptotic expansion of the fracture aperture, w That is why even very good coincidence of the 
crack widths cannot guarantee similar effect for the net fluid pressure. This also emphasizes importance 
of preserving further asymptotic terms of the solution for the quality of computations. 

The relative deviations of the approximations given in m from our numerical results are presented 
in Fig. 23 To complement this picture, we also present the case with n = 1 (Newtonian fluid) discussed 



Figure 22: KGD model (toughness dominated regime) - comparison of our numerical solution for N = 
100 with semi-analytical solution presented in (TT|: a) the crack opening w for n = 0, b) the crack 
opening w for n = 0.5, c) the fluid pressure p. 


It can be seen that we have very good coincidence of the analyzed results for n = 0. The graphs of 
w and p are indistinguishable from each other. However, this is no longer the case for n = 0.5, where 
the net fluid pressure curves diverge. The data for the case n = 1 is again very consistent with those 
obtained by the universal algorithm in [69]. 
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Figure 23: KGD model (toughness dominated regime) - comparison of the semi-analytical solution 
presented in El with our accurate numerical results for N = 100: a) relative deviation of the crack 
opening, w, b) absolute error of the fluid pressure, p. 


In order to explain this peculiarity, we should refer to the way in which the solution from E is 
reconstructed, and how the scaling between the respective equations in our formulation and that utilized 
in E is executed. In the recalled paper the author introduces a dimensionless crack half-length. This 
parameter is present in the leading terms of the power series representations for both the crack opening 
and the net fluid pressure, and thus in order to recreate the solution one needs to know its value. 
Moreover, the dimensionless crack half-length is employed in scaling when transforming the governing 
ODE from the formulation given in E to the formula ( |67| . In fact, the scaling criterion relates this 
value to the crack propagation speed 0(1). 

In general, there are two alternative methods of obtaining the dimensionless crack half-length: one 
utilizing the integral of w over the spatial interval (the fracture volume), and the other based on the 
fluid influx. Additionally, for n = 0 and n = 1 explicit formulae relating this parameter to one of the 
coefficients given in E can be found. For an exact solution, all methods should produce the same 
results, but naturally one cannot expect such an effect when using an approximation. Unfortunately, 
the author does not provide the numerical values of the dimensionless crack half-length. In [691 . when 
analyzing the data for a Newtonian fluid, the aforementioned explicit formula was used in computations. 
The same approach was implemented for n = 0. It is evident that such a method produces far better 
coincidence of the tip asymptotes between our results and the data from E. which leads to appreciably 
better convergence between the corresponding solutions than in the generic case (0 < n <1). Applying 
the alternative methods resulted in reduced solution accuracy for both, n = 0 and n = 1. For n = 0.5 we 
computed the dimensionless crack half-length from the non-local relation (based on the fracture volume). 
We believe that it constitutes a more reasonable choice than the influx criterion, as the global condition 
is less sensitive to the point-wise deterioration of the accuracy. It is noteworthy that, when computing 
the sought for parameter from the influx, its value differs from the previous one by approximately ten 
percent. 


4 Solution of the problem in transient regime 

In this section we discuss an extension of the algorithm presented above to the time-dependent variant 
of the problem. The main assumptions and blocks of the algorithm remain the same. The new fea- 
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tures introduced here are subroutines for approximating the temporal derivative and the crack length. 
Presented numerical examples demonstrate the performance of the general variant of the algorithm. 


4.1 Problem formulation 


Analogously to that done for the self-similar variant of the problem, let us write the basic system of 
equations. 

We use the following representation of the temporal derivative of the crack opening, as was proposed 
in [551: 


dw 

~dt 


= 2 


w(U+i, x) — w(U, x) dw 


ti +1 ti 


dt 


(94) 


The basic advantage of the presented formula over the ordinary finite difference is that it provides higher 
order approximation 0(At 2 )) than the latter (O(At)). Note that the value of the derivative at the initial 
time step can be computed analytically from the continuity equation (see details in [55]). 

The basic set of equations employed in the universal algorithm for the transient regime is listed below: 


• equation (941 to compute the temporal derivative, 


equation defining the reduced velocity (511 in the following transformed version: 



dw ( C A C(w)) 1/n \ 

L m,n + 1 ” + 


(95) 


solvability condition resulting from (95) and (52) utilized to compute the parameter C(w)\ 

f 1 dw(t,x) , qo(t) (C A £(w)) 1/,n f 1 f 1 

/ -- ~dx - —J + - 7 — ,, , , / w(t,x)dx+ / qi(t,x)dx~ 0, 

Jo dt L(t) L"*/n+i(t) J 0 v ; Jo 


(96) 


• equation to compute the crack opening, obtained by combining (55) with the respective inverse 
elasticity operator (26) (28): 


w = B 


L(t) <£( 77 ) + 1 ) 


{C A C{w)) l,n ' 

L m / n [t) 


(97) 


the influx boundary condition (52) and the regularized conditions at the crack tip, 


• initial conditions, which also allow computation of the temporal derivative at the initial time step, 


relation to compute the crack length (45). 


4.2 Computational algorithm for the transient regime 

The solution to the transient variant of the problem is computed at each time step by the iterative 
algorithm and is an extension of that presented for the self-similar formulation. By analogy to the 


description given in Section 3.3 we can define the following stages of computations at time ti+ 1 : 


• Preconditioning. The first approximation of the crack opening ^ (t i+1; x) is specified based on 
the temporal derivative from the previous time step w' t (ti,x ) and the initial condition w(ti,x). 
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• First step. The temporal derivative of the crack opening is computed using ([94]). Note that when 
obtaining the final solution w(ti+ 1 , x), one automatically has its temporal derivative as well. Then, 
we calculate the crack length L^~ 1 \ti + 1 ) using formula ( |45| with x)) (or C(w(U, x)) 

at the first iteration). Next equation ( f96j ) yields ( w(t i+1 , x)) which, substituted into fl95| , serves 
to compute the reduced particle velocity (i; + i, x). The e-regularization technique is applied to 
carry out the integral in (95). As a result, functions C < -^(w(U + 1 , z)) and ^ 3 \t i+ 1 , x) computed 
at this stage satisfy: i) fluid balance equation ( |96| ), ii) continuity equation ( |95| , iii) regularized 
boundary condition for <fi which is an equivalent of ( |53[ ), iv) the influx boundary condition (52) 
indirectly through the fluid balance equation. 


• Second step. Next the approximation of the crack opening (tj+i, x) is obtained from (97). The 
technique of numerically computing the operator B w is exactly the same as for the self-similar 
variant of the problem. 


The aforementioned two steps of the algorithm are repeated until respective components of the 
solution have converged within a prescribed tolerance. 

Remark 6. The modular algorithm architecture enables one to easily introduce the subroutine for 
crack length computation as an additional block. Naturally, this block was not present in the self-similar 
variant of the algorithm. 


4.3 Algorithm performance in the transient regime 


In this part of the paper we present a brief investigation of the algorithm performance in the transient 
regime. The few examples shown below constitute a demonstration of the potential of the proposed 
numerical scheme. This analysis is not exhaustive, as such a task would massively increase the volume 
of the paper (it would require e.g. presenting a number of simulations for various spatial meshing 
variants and time stepping strategies conducted for various values of the fluid behaviour index). As 
the influence of the spatial meshing on computational accuracy was thoroughly analyzed for the self¬ 
similar formulation of the problem, we will not consider this issue here. All results presented below were 
obtained for an adaptive spatial mesh composed of N = 100 nodal points, which for the self-similar 
variant produced a maximal solution error of at most order 10~ 6 for both the crack opening and particle 
velocity, regardless of the HF model. It is obvious that in the transient regime one cannot expect the 
same level of accuracy. Indeed, as was shown in |69j . the solution accuracy is naturally determined by 
the relation between the spatial and temporal meshing. For each of the HF models we present the results 
of computations for two values of the fluid behaviour index (n = 0.3 and n = 0.7) and two different 
temporal meshes. The number of time steps, denoted by M, was set to M = 20 and M = 100, while 
the same normalized target time t = 100 was assumed for all computations. The time meshing strategy 
which provides an approximately parabolic time step distribution (with finer meshing for small times) 
was adopted from [68] . 

The benchmark solutions utilized here are the time-dependent extensions of self-similar benchmarks 


used in Sections 3.4.1 (PKN), 3.4.2 (viscosity dominated regime of KGD) and 3.4.3 (toughness dominated 
regime of KGD). The analyzed error measures are: the relative error of crack opening, Sw(t,x), the 
relative error of the particle velocity, 5v(t,x ), and the relative error of the crack length, SL(t). 


4.3.1 Analysis of the algorithm - PKN model 

In Fig. [24] - [25] the distributions of relative errors over space and time for the crack opening, 6w, and 
the particle velocity, Sv, are shown for n = 0.3. Regardless of the temporal mesh density we see some 
general trends. Namely, for a fixed number of the time steps, M, the graphs for 5w and Sv are practically 
identical. Almost uniform distribution of the computational errors over the crack length, x, is observed. 
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When analyzing the time characteristic it is evident that the greatest error growth takes place for small 
times before the level of relative error stabilizes. In the case under consideration, by taking M = 100 
instead of M = 20, one reduces the error by over an order of magnitude. One should remember that 
M = 20 produces a very coarse temporal mesh (only 20 steps from the initiation of fracture growth to 
the large time asymptote). 
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Figure 24: PKN model. Relative solution error for n = 0.3 with N = 100 (spatial mesh), M = 20 
(temporal mesh): a) the crack opening Sw, b) the particle velocity Sv. 



Figure 25: PKN model. Relative solution error for n = 0.3 with N = 100 (spatial mesh), M = 100 
(temporal mesh): a) the crack opening Sw, b) the particle velocity Sv. 

A bit different trend is observed in Fig. [26] [27j where the relative errors for n = 0.7 are presented. 
First, the overall accuracy has improved by approximately one order of magnitude with respect to the 
case of n = 0.3. The maxima of Sw are now about two times smaller than the maxima of Sv for the 
same M. Also, the potential for error reduction when increasing M is larger. 

When analyzing the crack aperture it shows that the transition from M = 20 to M = 100 gives two 
orders of magnitude better accuracy. In the case of particle velocity this reduction is smaller, but still 
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distinctly greater than it was for n = 0.3. Nonetheless, it is difficult to speculate whether the reason 
for this change in general trends can be entirely attributed to changing the value of n. Note that when 
expanding the self-similar benchmarks to their time-dependent versions they produce quantitatively 
different dynamic behaviour (temporal derivatives, crack propagation speed etc.). 



Figure 26: PKN model. Relative solution error for n = 0.7 with N = 100 (spatial mesh), M = 20 
(temporal mesh): a) the crack opening Sw, b) the particle velocity Sv. 



Figure 27: PKN model. Relative solution error for n = 0.7 with N = 100 (spatial mesh), M = 100 
(temporal mesh): a) the crack opening Sw, b) the particle velocity Sv. 


To complement this part of the analysis let us look at Fig. 36i, where the crack length error, SL, 
for the respective cases is depicted. It shows that, regardless of the value of fluid behaviour index n and 
the number of time steps M , the error level stabilizes after the initial growth, however its magnitude 
depends on the considered variant. In general, a better solution for w and v produces better results for 
L , which is an intuitive trend. 
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4.3.2 Analysis of the algorithm - KGD model in viscosity dominated regime 


Let us now analyze the accuracy of computations for the viscosity dominated regime of the KGD model. 
In Fig. 28 - 29 the relative errors of w(t,x) and v(t, x) for n = 0.3 are presented. This time (unlike in 


the PKN variant) the maximal values of Sv are approximately two times greater than the maxima of 
Sw. Moreover, the distribution of Sw exhibits different shapes than that of Sv. Increasing the number of 
time steps from 20 to 100 entails a two order improvement in accuracy for the crack opening, and one 
order for the particle velocity. In comparison with the PKN case (Fig. 24 - Fig. 25), the level of Sw 
increased approximately two times, while the magnification of Sv was even greater. The latter effect is 
caused by the error increase at the crack tip for the small time range. 



Figure 28: KGD model (viscosity dominated regime). Relative solution error for n = 0.3 with N = 100 
(spatial mesh), M = 20 (temporal mesh): a) the crack opening Sw, b) the particle velocity Sv. 



Figure 29: KGD model (viscosity dominated regime). Relative solution error for n = 0.3 with N = 100 
(spatial mesh), M = 100 (temporal mesh): a) the crack opening Sw, b) the particle velocity Sv. 

Fig. [30] [TT| depict the solution error for n = 0.7. This time, for fixed M, we obtain the same level 
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of inaccuracy for both the crack opening and the particle velocity. In this case the magnification of Sv 
at the crack tip is not observed. Both errors remain stable over the whole range of x. In comparison 
with the case when n = 0.3, the accuracy is of at least one order of magnitude better. 



Figure 30: KGD model (viscosity dominated regime). Relative solution error for n = 0.7 with N = 100 
(spatial mesh), M = 20 (temporal mesh): a) the crack opening Sw, b) the particle velocity Sv. 
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Figure 31: KGD model (viscosity dominated regime). Relative solution error for n = 0.7 with N = 100 
(spatial mesh), M = 100 (temporal mesh): a) the crack opening Sw, b) the particle velocity Sv. 


The relative errors for the crack length L{t) are presented in Fig. 36 3. The character of growth for 
SL, as well as its magnitude, are very similar to those reported previously for the PKN model (Fig. 36 1 ). 
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4.3.3 Analysis of the algorithm - KGD model in toughness dominated regime 


The results of computations for n = 0.3, illustrated by the relative solution error, are presented in 
Fig. [32] [33] The level of Sw for different M is now one order of magnitude lower than that for the 
viscosity dominated benchmark, which we believe is not a feature of the algorithm itself but rather of 
the analyzed benchmark (as was shown in the self-similar formulation, the solution accuracy depends 
on the benchmark type). In the case of particle velocity, the magnitude of maximal errors remains 
comparable to those obtained previously in the viscosity dominated regime, but the spacing of Sv is 
completely different. This time a distinct error magnification is observed at the crack tip over the whole 
time interval. It can be explained by the fact that the crack propagation speed is now computed based on 
the two leading asymptotic terms of w (instead of one for the PKN model and KGD model in viscosity 
dominated regime, compare (42) - ©). The multiplier of the second term is approximated, by its 
very nature, in a less accurate way than the coefficient of the leading term, which may affect the overall 
accuracy of v(t, 1). 
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Figure 32: KGD model (toughness dominated regime). Relative solution error for n = 0.3 with N = 100 
(spatial mesh), M = 20 (temporal mesh): a) the crack opening Sw, b) the particle velocity 5v. 

The results for n = 0.7 are shown in Fig. [34]- Fig. [35] The graphs for Sw resemble those obtained for 
n = 0.3 in both the level and distribution of errors. By taking M = 100 instead of M = 20 one increases 
the accuracy by two orders of magnitude. In the case of particle velocity we do not have a characteristic 
error increase in the near-tip region, which was present for n = 0.3. However, a sharp magnification of 
Sv(t, 1) for small time is still in place. Also, some minor fluctuations in the tip velocity error for larger 
times are observed. Compared to the variant with n = 0.3 the maximal level of Sv remained the same 
for M = 20, but was reduced by one order of magnitude for the refined temporal meshing (M = 100). 

Finally, the graphs for the crack length error are shown in Fig. |36( c. This time the differences in the 
value of SL are much smaller than that for other HF models. Moreover, unlike the previous cases, now 
it is n = 0.3 which gives slightly better results. However the general trend of error increase in the early 
time, and its stabilization after the initial growth, is still present. 


5 Discussion and Conclusions 

In the paper a problem of ID hydraulic fracture has been analyzed. The power law rheological model 
of the fracturing fluid has been employed for the classical PKN and KGD (in viscosity and toughness 
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Figure 33: KGD model (toughness dominated regime). Relative solution error for n = 0.3 with N = 100 
(spatial mesh), M = 100 (temporal mesh): a) the crack opening Sw, b) the particle velocity Sv. 
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Figure 34: KGD model (toughness dominated regime). Relative solution error for n = 0.7 with N = 100 
(spatial mesh), M = 20 (temporal mesh): a) the crack opening Sw, b) the particle velocity Sv. 


dominated regimes) models. Following the idea from [69| the problem has been reformulated in terms 
of a new pair of dependent variables: the crack opening, re, and the reduced particle velocity, cj). A 
unified equation (411 defining the interrelation between the crack propagation speed and the solution tip 
asymptotics has been derived for various elasticity operators and crack propagation regimes. Self-similar 
formulations of the problem for different HF models have been delivered. 

A universal algorithm for solving the problem, being an extension of that proposed in [69] , has 
been introduced. It enables us to tackle various elasticity operators, fluid flow and fracture propagation 
regimes in the framework of a unified numerical scheme. The computational algorithm is comprised of 
two basic modules: 


• Universal module to compute the reduced particle velocity - this block remains the same regardless 
of the HF model and crack propagation regime. 
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Figure 35: KGD model (toughness dominated regime). Relative solution error for n = 0.7 with N = 100 
(spatial mesh), M = 100 (temporal mesh): a) the crack opening 6w, b) the particle velocity Sv. 
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Figure 36: Relative error of the fracture length: a) PKN model, b) KGD model in viscosity dominated 
regime , c) KGD model in toughness dominated regime. Spatial mesh: N = 100. 


• Module to compute the crack opening - this block is adjusted depending on the elasticity relation 
in use. The asymptotic parameters embedded here facilitate adaptation of the module to various 
crack propagation regimes. 


Crucial elements of the proposed algorithm are: 

• The mechanism of fracture front tracing based on the local condition of Stefan type (speed equa¬ 
tion). This device utilizes an explicit formula (411 relating the crack propagation speed to the 
crack opening tip asymptotics. This in turn allows us to compute in an explicit form the crack 
length from the equation (45). 

• Appropriate handling of the independent (spatial and temporal meshing) and dependent variables. 

• Improved approximation of the temporal derivative. 

• Dedicated regularization techniques (the so-called e-regularization and operator regularization of 
the governing equations). 
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The algorithm performance has been thoroughly investigated. To this end a set of analytical bench¬ 
marks, introduced in Appendix A, have been utilized. The accuracy of computations has also been 
verified against other solutions available in the literature. The following conclusions can be drawn about 
the performance of the proposed numerical scheme: 


• The algorithm is efficient and stable regardless of the hydraulic fracture model and crack propa¬ 
gation regime. 

• High solution accuracy can be obtained at low computational cost. Proper utilization of the basic 
devices embedded in the algorithm allows good results to be obtained for even coarse meshing. 

• Even in the extreme regimes of crack propagation (e.g the small toughness case) high solution 
quality can be retained, however it necessitates appropriate measures such as skilled handling 
of the e-regularization technique and spatial meshing. As is to be expected, the computational 
efficiency decreases for such limiting variants of the problem. 

• The numerical results provided by the computational scheme are much more accurate than any 
other data available in the literature (e.g. for 100 spatial mesh points at least three orders better 
accuracy can be obtained regardless of the analyzed HF model). 


The developed numerical algorithm also allowed us to improve two important classical benchmark 
solutions. For the PKN model the formula (87) for the self-similar crack half-length, £*, when combined 
with the solution given in m provides an improved semi-analytical benchmark for an arbitrary value 
of the fluid behaviour index. Additionally, a new improved solution approximation of the viscosity 
dominated regime of the KGD model has been proposed. The approximation guarantees much better 
accuracy than any other semi-analytical solution available in the literature. Its quality is equally good 
for the whole class of shear-thinning fluids. The simple form of the proposed formulae facilitates their 
application in testing other numerical algorithms. 

To summarize, the proposed algorithm has been proven to posses all of the advantages of the basic 
version originally introduced in [69] for the Newtonian fluids. It constitutes a flexible and reliable tool 
for numerically modeling hydraulic fractures driven by shear-thinning fluids. 
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Appendix A: Benchmark solutions 


Let us describe a set of benchmark solutions used to analyze the accuracy of computations. The main 
concept for their construction was already presented and employed in [55]. The basic assumption is 
that, if we can provide an analytical solution for the self-similar problem defined in Section |3j it can be 
immediately extended to the time-dependent form by relations (56). 

The general idea of finding the benchmark solution is as follows. First, let us assume that the crack 
opening function can be expressed as a weighted sum of appropriately chosen base functions: 


R 

w(x) = ^2 A iClifx). 

i =0 


(98) 
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The functions 12; (x) are selected in a way which enables one to: i) comply with the respective asymp¬ 
totic representation (38), ii) compute the pressure operators (58) analytically, iii) satisfy the boundary 
conditions (69) (and (70) for KGD). Provided that ii) is fulfilled, the pressure function can be calculated 

R 


in a closed form from (58): 


P(x) =^2Xiili(x), 


(99) 


i =0 


where each of the functions II; (x) corresponds to respective function 12; (x). Then, according to (63), 
the particle velocity function gives: 


v{x) = 


n+1 


R 


- ^AA(x) £> 


; i—0 


i=0 


dUj 

dx 


( 100 ) 


The reduced velocity can then be determined by using (100) in (65): 

\ n+1 


4>(%) = 


( R \ 1 ' x R 

^A;12;(x)J J^A;^ - x(c A C(w)^ ". 

2=0 / 2=0 


( 101 ) 


Next, by applying (98) and (101) in (67), the benchmark leak-off function qi(x) can be defined. The 
benchmark value of qo is determined by substituting (98) and (101) into the boundary condition ( 69 ) 3 : 


do = 


R 


n+2 


R 


diii , 




V. 2=0 


2=0 


( 102 ) 


In this way, the analytical benchmark solution is fully defined by (98), (101), the corresponding leak-off 
function and the corresponding influx value (102). The respective fluid balance equation <zu is satisfied 
automatically. 

Al: PKN model 

For the PKN model let us adopt the following base functions rc;(x): 

£li{x) = (1 -+++S i = 0,..,R— 1, n R {x) = e x (l -x) fl+ +3. (103) 

The function fl R (x) was taken in this manner in order to introduce an additional non-linear effect to 
the benchmark, without violating the asymptotic behaviour. 

By applying representation (103) in (98) and elasticity operator (59), one obtains the formula for the 
pressure function, p(x), which after differentiation yields: 

, rR-i 


= ~ 


X> 

,i=0 


1 


n + 2 


(1 - xf + X R e x (l - x) R+ s+5 


R 


1 — x 


- 1 


(104) 


Then by formulae (100) (102) one can construct the benchmark solution, taking C{w) = Ag +2 . 
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A2: KGD model 

For this model we assume that the self-similar crack aperture is defined by the following base functions: 

1 + \/l — x 2 


fli(x) = (l-x 2 )“ i C'“ ( t + 1 1 ^ 2 (a;), i = 0, fi R (x) = \/l - x 2 --(l-x 2 ) 3/2 -x 2 In 

(105) 

where y_ 2 (x) is the ultraspherical polynomial. Respective powers at are to be taken in accordance 
with the asymptotic behaviour of the crack opening (see e.g. Table [l]). The term f 1r(x) was introduced 
to obtain a non-zero pressure gradient for x = 0. Note that: 


Cl R (x) = O ((1 


-- 2 ) 5/2 ), 


1 . 


(106) 


By applying (105) in (58)2 and subsequent differentiation one obtains an analytical formula for the 

R -1 


pressure gradient: 


A 

dx 


P( x ) = - X] X i H i( X ) - ^R - X ) > 


2 — 0 


where: 


H z (a.i,x) = 

3 

2 Fi(- - i - a 


H 0 (x) = -a 0 (2a 0 - l)B(l/2, a 0 )x - 2 Fi(3/2 - a 0 , 2; 3/2; a: 2 ), 


(1 — 2a*)(2i + 1) 


xB( — - — i,a + i - |-l)-| [2(i + 1 )(* + 1 + a) — 2 — 5i] x 
,* + 1; A,x 2 ) - 2 ^{i + l)(2i -3 + 2 a)x 2 ■ 2 +i(^ -i-a,i + 2; * 


47 t (* + a) 
3 


(107) 

(108) 
(109) 


> 1 . 


Here B is the beta function and 2 +l is the Gauss hypergeometric function. Then from ( |100 ) 
( |102 ) we have the complete benchmark solution, where C(w) = Ag +2 (viscosity dominated regime) 
or C(w) = Aq +1 Ai (toughness dominated regime). Moreover, by properly handling the values of re¬ 
spective coefficients Aj, one can control the asymptotic behaviour of the particle velocity and leak-off 
function, and in the case of the toughness dominated mode the value of the self-similar stress intensity 
factor Kj. 


Remark 7. For the toughness dominated regime of the KGD model we take: 

Do(^) = \/l — x 2 . 


( 110 ) 


As a result the corresponding pressure component is constant along the entire fracture length which 
implies that: 

H 0 {x) = 0. (Ill) 

Remark 8. If need be one can also utilize a simplified version of the benchmark solution. Namely, 
when functions £li(x) from the basis (105) are taken as the power terms fli(x) = (1 — x 2 ) ai , respective 
entries Hi(x) from (107) are then computed according to (108). All other benchmark components are 
determined in the same way as previously. 

A2.1: Perfectly plastic fluid (n=0) 


In the case of a perfectly plastic fluid (n = 0) the Poiseuille equation degenerates to the form (78). 
For this reason it can no longer be used to compute the particle velocity, and as such the method of 
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constructing benchmark examples described above fails. In order to circumvent this deficiency we shall 
modify the problem slightly and propose a set of analytical benchmark solutions for such a formulation. 
As before we take some representation of the crack opening (|98[) which accounts for the proper tip 


asymptotics. When substituting it into the elasticity relation (|58|)2 one arrives at the expression for the 


net fluid pressure in the form (99). Next we impose the following constraints on the values of the first 
coefficients A,-: 


• the viscosity dominated regime of the KGD model 

Ao = 

• the toughness dominated regime of the KGD model 

4 

AoAl = -^- 


As a result of these assumptions the equivalent of (78) yields: 

dp 


-w—=g(x), g(l) = 1 , 
ax 


( 112 ) 


(113) 


(114) 


where g(x) is a known function depending on the chosen w and p. 


In order to construct the set of benchmark solutions we assume that (1141 holds instead of (78). 


Consequently, the computational formula (81) needs to be modified to: 


w{ x) = — 


l ^r\ K ( r i’ x )d'n + -^=.Ki-v/i -: 
Jo W(V) V 71 " 


(115) 


In this way, by making slight modifications to the problem (introduction of relations (114) and ( 115| ) 
instead of (781 and ©) we arrive at the formulation which allows one to find a class of analytical 
benchmark solutions that can be used to test the numerical algorithm. 


The benchmark value of the self-similar crack propagation speed Do is computed according to (82), 


with the leak-off function and influx magnitude chosen for convenience. Finally, the reduced particle 


velocity is obtained in an analogical way from (73). 


Throughout this paper we use the following variants of the benchmark solutions for n = 0. 

• The viscosity dominated regime. 

Three base functions fli are employed: 

Cl 0 (x) = 1 — x 2 , Cli(x) = (1 — x 2 ) 2 , f) 2 (x) = CIr(x), 

where f Ir is defined in (|105). The corresponding components of the pressure derivative are com¬ 


puted according to ( 1071) (108). 


• The toughness dominated regime. 

The base functions f \ are taken as: 

&o(x) = \J\ — x 2 , Di(x) = (1 — x 2 ) 3 / 2 log(l — x 2 ), 


D 2 (x) = (Ir(x), 


where Qr is defined in (105). As usual, the square root term produces a constant pressure com¬ 


ponent (zero gradient). The remaining components of the pressure derivative are: 

,, . . ( 5 A 3 1 — 2x 2 . 7r 

Pi(x) = (^3log2 - -J x - - arcsm(x) —j ===, p 2 (x) = x - 

Note that for the foregoing representation of the crack opening, the integrals needed to define the 
benchmark values of particle velocity and reduced particle velocity are computed analytically. 
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Appendix B: Coefficients of the formulas (90) — (92) 


Coefficients wi(n), pi(n) and vi(n) in formulas (90) (92) are all functions of the fluid behaviour index 


n and are approximated by the following expression: 


tf(rc) = X/ 




(116) 


/c—0 


with values of rf given in Table |4j 


i?(n) 

r 0 

rf 

r # 

'2 

'3 

r # 

7 4 

Wi 

-0.165 

0.529 

-0.7587 

0.6006 

-0.1985 

W 2 

0.673 

-0.4886 

0.00747 

0.1486 

-0.0692 

Pi 

0.6228 

-0.1334 

-0.5781 

0.5167 

-0.1536 

P2 

-1.388 

0.1244 

1.977 

-1.873 

0.5961 

P3 

1.911 

-1.521 

-0.4046 

0.8542 

-0.2925 

PA 

-1.129 

1.952 

-1.674 

0.9952 

-0.3012 

P5 

0.04245 

-0.1925 

0.2071 

-0.1114 

0.02931 

P6 

0.05075 

-0.1895 

0.2868 

-0.2337 

0.07845 

V 1 

0.00987 

-0.1373 

0.1017 

-0.0497 

0.0122 

V 2 

0.01939 

-0.00148 

0.02551 

-0.01331 

0.00323 

h 

0.2264 

-0.06635 

-0.04062 

0.03081 

-0.00845 

VA 

-0.3275 

0.1403 

0.04508 

-0.04654 

0.01473 

i>5 

0.2098 

-0.09169 

-0.01172 

0.01446 

-0.00373 


Table 4: The values of coefficients rt in equation (116). 
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Nomenclature 


Kj stress intensity factor 

Kic material toughness 

l{t) crack length (L(t) after normalization) 

n fluid behaviour index 

p(t , x) net fluid pressure 

q{t , x) fluid flow rate 

qi(t,x) leak-off 

qo{t) influx 

v(t,x) particle velocity 

w(t,x) crack opening 

qzq , «i, 02 exponents of leading terms in asymptotic expansion of the crack opening 

Pi, P 2 exponents of leading terms in asymptotic expansion of the particle velocity 

Co> Ci; C2 exponents of leading terms in asymptotic expansion of the reduced particle velocity 

<j>(t, x) reduced particle velocity 

All variables with ” symbol refer to the normalized formulation, while ” A” symbolizes self-similar 
formulation. 


53 



